The hydraulic efficiency of single fractures: correcting the cubic law parameterization for self-affine surface roughness and fracture closure

Abstract

Quantifying the hydraulic properties of single fractures is a fundamental requirement to understand fluid flow in fractured reservoirs. For an ideal planar fracture, the effective flow is proportional to the cube of the fracture aperture. In contrast, real fractures are rarely planar, and correcting the cubic law in terms of fracture roughness has therefore been a subject of numerous studies in the past. Several empirical relationships between hydraulic and mechanical aperture have been proposed based on statistical variations of the aperture field. However, often, they exhibit non-unique solutions, attributed to the geometrical variety of naturally occurring fractures.

In this study, a non-dimensional fracture roughness quantification scheme is acquired, opposing effective surface area against relative fracture closure. This is used to capture deviations from the cubic law as a function of quantified fracture roughness, here termed hydraulic efficiencies. For that, we combine existing methods to generate synthetic 3-D fracture voxel models. Each fracture consists of two random, 25 cm${}^{\mathrm{2}}$ wide self-affine surfaces with prescribed roughness amplitude, scaling exponent, and correlation length, which are separated by varying distances to form fracture configurations that are broadly spread in the newly formed two-parameter space (mean apertures in submillimeter range).

First, we performed a percolation analysis on 600 000 synthetic fractures to narrow down the parameter space on which to conduct fluid flow simulations. This revealed that the fractional amount of contact and the percolation probability solely depend on the relative fracture closure.

Next, Stokes flow calculations are performed, using a 3-D finite differences code on 6400 fracture models to compute directional permeabilities. The deviations from the cubic law prediction and their statistical variability for equal roughness configurations were quantified. The resulting 2-D solution fields reveal decreasing cubic law accordance down to 1 % for extreme roughness configurations. We show that the non-uniqueness of the results significantly reduces if the correlation length of the aperture field is much smaller than the spatial extent of the fracture. An equation was provided that predicts the average behavior of hydraulic efficiencies and respective fracture permeabilities as a function of their statistical properties. A model to capture fluctuations around that average behavior with respect to their correlation lengths has been proposed. Numerical inaccuracies were quantified with a resolution test, revealing an error of 7 %.

By this, we propose a revised parameterization for the permeability of rough single fractures, which takes numerical inaccuracies of the flow calculations into account. We show that this approach is more accurate compared to existing formulations. It can be employed to estimate the permeability of fractures if a measure of fracture roughness is available, and it can readily be incorporated in discrete fracture network modeling approaches.

1 Introduction

The geometrical inhomogeneities of single fractures and their effect on fluid flow remain a crucial parameter for understanding the hydraulic properties of fractured reservoirs, such as crystalline or tight carbonate rocks with nearly impermeable matrices. Hence, it has wide-ranging industrial applicability in the fields of petroleum and gas production, geothermal energy recovery, ${\mathrm{CO}}_{\mathrm{2}}$ sequestration, nuclear waste disposal, and groundwater management. Fluid flow in fractured reservoirs is commonly modeled by the discrete-fracture-network (DFN) approach , which relies on knowing the permeability of single fractures. The permeability of a single fracture is often approximated by the well-known cubic law , assuming that a fracture is composed of two parallel plates separated by a constant aperture. However, natural fracture walls show deviations from planarity, i.e., roughness, resulting in varying apertures within the fracture plane. On top of that, fluid–rock interactions like dissolution , erosion , and mineral growth , as well as the surrounding stress field further modify the geometry of a fracture, causing deviations of the parallel plate assumption.

Considerable effort has been made to study the effect of fracture surface roughness on flow and reactive transport behavior. Early attempts employed the 2-D Reynolds equation, a simplification of the Navier–Stokes equations, which assumes that the cubic law holds locally with the aperture varying in the $x$–$y$ along-fracture plane. They derived semi-empirical functions that describe the deviations from the cubic law in terms of the mean and standard deviation of the aperture field. Increasing computational power led to numerical improvements, with 3-D lattice Boltzmann or Navier–Stokes simulations revealing the non-uniqueness of previous functional approximations of fracture permeability. Factors such as shear displacement , tortuosity, and the degree of mismatch between the opposing fracture surfaces were demonstrated to affect fluid flow paths and permeabilities. Detailed analyses of exposed fracture surfaces have shown that self-affine fractal models are capable of quantifying surface roughness properties from thin section scale to outcrop scale. Thereby, the dependence of surface roughness as a function of observation scale is captured by their scaling exponent (the so-called Hurst exponent, $H$). For example, mode I fractures in a porous sandstone showed $H$ of 0.4–0.5 . Microfractures in Pomeranian shale featured $H$ of $\mathrm{0.3}$ and $\mathrm{0.5}$ , depending on the opening mode. In other studies, $H$ of fault surfaces (mode II fractures) tends to fall in the range of 0.6–0.8 with respect to slip orientation , regardless of rock type. Based on this, it is commonly assumed that a fracture consists of two opposing self-affine surfaces, and the resulting aperture field follows the same scaling relationship, assuming both surfaces are uncorrelated . However, observations of opposing fracture walls have demonstrated that the two surfaces tend to be well correlated above a specific length scale and non-correlated below it, which poses an upper limit to the self-affine scaling in nature. Following , the ratio between system size ($L$) and the correlation length ($l_{\mathrm{c}}$) defines whether the fracture has an intrinsic permeability or not. Their statistical approach suggested that permeabilities of uncorrelated fractures (i.e., $l_{\mathrm{c}}/L=\mathrm{1}$) are strongly fluctuating and anisotropic for the same roughness configurations, revealing the importance of considering low $l_{\mathrm{c}}/L$ ratios to be able to quantify an intrinsic fracture permeability.

Although extensively studied, no clear mathematical relationship between fracture roughness and permeability has been derived so far, leaving the cubic law as standardized parameterization in DFN modeling approaches. Thus, an applicable refinement is desired to promote their realism to help better understanding fluid flow on a reservoir scale. In this paper, existing algorithms are used to generate a large dataset of synthetic fractures covering all possible kinds of roughness configurations. Single-phase 3-D Stokes flow calculations are then performed with a finite difference code, utilizing a high-performance-computing (HPC) cluster to handle the associated computational effort. By interpreting the statistical variations of the results, a refined parameterization of single-fracture permeability is proposed, which is demonstrated to provide accurate predictions for the permeability of rough fractures.

2 Method and data

2.1 Fluid flow in self-affine fractures

Generally, the flow of an incompressible Newtonian fluid is most accurately described by the Navier–Stokes equations (NSEs). Assuming that the flow is solely laminar Reynolds numbers below unity according to, the fluid viscosity is constant, and gravity is negligible at the system size, they reduce to the simpler Stokes equations, i.e., momentum balance (Eq. ) and continuity (Eq. ) equations, which for steady-state flow conditions are given in compact form by $$\begin{array}{}\text{1}& {\displaystyle }\mathit{\mu }{\mathrm{\nabla }}^{\mathrm{2}}\mathit{v}=\mathrm{\nabla }P,\text{2}& {\displaystyle }\mathrm{\nabla }\cdot \mathit{v}=\mathrm{0},\end{array}$$ with the fluid's dynamic viscosity $\mathit{\mu }$, pressure $P$, and velocity vector $\mathit{v}=(v_x,v_y,v_z)$, $\mathrm{\nabla }$, $\mathrm{\nabla }\cdot$, and ${\mathrm{\nabla }}^{\mathrm{2}}$ denote the gradient, divergence, and Laplace operator for 3-D Cartesian coordinates, respectively.

These equations can be solved analytically for an idealized fracture, consisting of two parallel plates, vertically separated by a constant aperture $a$. Volumetric integration of the horizontally extended Poiseuille-flow solution yields the well-known cubic law:

3$$Q=-{\displaystyle \frac{w{a}^{\mathrm{3}}\mathrm{\Delta }P}{\mathrm{12}\mathit{\mu }}},$$ with total volumetric flow rate $Q$, fracture width $w$, and pressure gradient along the fracture $\mathrm{\Delta }P$ (see for a more detailed derivation). Combining Eq. () with Darcy's law for flow through porous media, 4$$Q=-{\displaystyle \frac{kA\mathrm{\Delta }P}{\mathit{\mu }}},$$ with cross-sectional area $A$, leaves the intrinsic permeability $k$ of an idealized fracture proportional to its aperture by $k\propto a^{\mathrm{2}}/\mathrm{12}$. For a rough walled fracture, the aperture is no longer constant but rather varying across the fracture plane. The mean planes of an upper and lower rough surface $s_{\mathrm{u}}(x,y)$ and $s_{\mathrm{l}}(x,y)$ are separated by a constant distance $d$ to form $a(x,y)$ according to 5$$a(x,y)=\left[s_{\mathrm{u}}(x,y)+{\displaystyle \frac{d}{\mathrm{2}}}\right]-\left[s_{\mathrm{l}}(x,y)-{\displaystyle \frac{d}{\mathrm{2}}}\right].$$ Depending on $d$, the surfaces may overlap at certain points and form contact areas, such that $a_{\mathrm{0}}(x,y)$ is zero at these locations: 6$$a_{\mathrm{0}}(x,y)=\left\{\begin{array}{ll}a& \text{if }a(x,y)>\mathrm{0}\\ \mathrm{0}& \text{if }a(x,y)\le \mathrm{0}.\end{array}\right.$$ Assuming that $s_{\mathrm{u}}$ and $s_{\mathrm{l}}$ are self-affine, the standard deviation $\mathit{\sigma }$ of the aperture field computed at system size $l$ scales as : 7$$\mathit{\sigma }=\left\{\begin{array}{ll}\mathit{\beta }{l}^H& \text{if }\mathrm{0}\le l\le l_{\mathrm{c}}\\ \mathit{\beta }{l}_{\mathrm{c}}^H& \text{if }{l}_{\mathrm{c}}\le l\le L,\end{array}\right.$$ with the maximal system size $L$ and the correlation length $l_{\mathrm{c}}$. Below $l_{\mathrm{c}}$, the fracture is uncorrelated, and it is well correlated above it (see Fig.  for a visual explanation). The prefactor $\mathit{\beta }$ delivers information about the overall amplitude of the surface roughness. $H$ typically denotes the scaling or Hurst exponent with $\mathrm{0}<H\le \mathrm{1}$, whereas $H=\mathrm{1}$ corresponds to self-similar and $H<\mathrm{1}$ to self-affine scaling e.g.,. Physically, self-affinity is expressed by higher height fluctuations at smaller scales, leaving $H$ as a measure for the ratio of large-scale versus small-scale roughness intensity.

Figure 1

(a) Aperture field of a synthetically generated fracture ($d=\mathrm{1}$ mm, ${\mathit{\sigma }}_{\mathrm{s}}=\mathrm{0.5}$ mm, $H=\mathrm{0.9}$) with a $l_{\mathrm{c}}/L$ ratio of $\mathrm{1}/\mathrm{16}$. (b) Zoomed-in image showing the uncorrelated part of the aperture field. Units of $L_x$ are in millimeters.

[Figure omitted. See PDF]

Figure 2

Two aperture fields constructed from synthetic fractures. Both aperture fields are based on the same sets of random numbers with varying Hurst exponents $H$, which are (a) $\mathrm{0.4}$ and (b) $\mathrm{0.8}$. The two statistical parameters ($\stackrel{\mathrm{‾}}{a}$ and ${\mathit{\sigma }}_{\mathrm{a}}$) are indicated by solid and dashed black lines, respectively. Axis units are in millimeters, while the vertical axis (indicating aperture) is exaggerated by a factor of 2 for clarity. Note that $\stackrel{\mathrm{‾}}{a}$ and ${\mathit{\sigma }}_{\mathrm{a}}$ are identical for panels (a) and (b). Increasing height fluctuations at smaller scales, caused by a lower Hurst exponent, results in a larger effective surface area $S$ for fracture (a) compared to panel (b).

[Figure omitted. See PDF]

Here, we use the following two non-dimensional quantities to quantify the geometry of a rough fracture: (i) the relative closure $R$ and (ii) the effective surface area $S$. We compute the relative closure by dividing the standard deviation of the aperture field at the maximal system size ${\mathit{\sigma }}_{\mathrm{a}}$ by the average aperture field $\stackrel{\mathrm{‾}}{a}$: 8$$R={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{a}}}{\stackrel{\mathrm{‾}}{a}}},$$ with $\stackrel{\mathrm{‾}}{a}$ defined by 9$$\stackrel{\mathrm{‾}}{a}={\displaystyle \frac{\mathrm{1}}{L_x{L}_y}}\underset{x=\mathrm{0}}{\overset{L_x}{\int }}\underset{y=\mathrm{0}}{\overset{L_y}{\int }}{a}_{\mathrm{0}}(x,y)\mathrm{d}x\mathrm{d}y.$$ This quantity or its reciprocal is commonly used to infer the amount of contact between the opposing fracture walls . Theoretically, it falls in the range $\mathrm{0}<R\le \mathrm{\infty }$, whereas $R=\mathrm{0}$ shows perfect accordance with parallel plates and the surfaces are in contact if $R\ge (\mathrm{3}\sqrt{\mathrm{2}}{)}^{-\mathrm{1}}$ see.

Furthermore, one requires a parameter that quantifies the effective surface roughness of a fracture since fractures with different $H$ can have equal $R$ values (see Fig.  for a visualization of the non-uniqueness of $R$). We therefore introduce a new quantity, the “effective surface area $S$”. This parameter uniquely combines varying amplitudes and scaling exponents, because an increase in fracture surface area is the direct consequence of increasing roughness. For that, we calculate the ratio of the surface area of the fracture ${\mathrm{sa}}_{\mathrm{f}}$ to twice the area of its projection on the fracture plane (i.e., 2 times the base area perpendicular to the flow direction) ${\mathrm{sa}}_{\mathrm{c}}$ and normalize it with the fractional amount of the aperture field that has opened, i.e., 10$$S=\left({\displaystyle \frac{{\mathrm{sa}}_{\mathrm{f}}}{{\mathrm{sa}}_{\mathrm{c}}}}\right)\left({\displaystyle \frac{\mathrm{1}}{\mathrm{1}-c}}\right),$$ with $c$ being the contact fraction of the aperture field (i.e., where $a_{\mathrm{0}}(x,y)=\mathrm{0}$), which leaves $\mathrm{1}\le S\le \mathrm{\infty }$, with $S=\mathrm{1}$, showing perfect accordance with parallel plates. To finally quantify the influence of fracture roughness on its intrinsic permeability, the proportionality resulting from the cubic law needs to be corrected e.g., by applying a correction factor according to 11$$k\propto \mathit{\chi }{\displaystyle \frac{a^{\mathrm{2}}}{\mathrm{12}}}.$$ The approximation of $\mathit{\chi }$ in terms of quantified fracture roughness, i.e., $\mathit{\chi }(R,S)$, is the main subject of this study.

2.2 Numerical permeability estimation

To simulate the laminar flow of an incompressible, isothermal, and isoviscous fluid, we use a 3-D binary voxel model input. We solve the linearized momentum balance (Eq. ) and continuity (Eq. ) equations in 3-D, using velocity and pressure as primary variables. The coupled system is implemented in the open-source software package LaMEM (Lithosphere and Mantle Evolution Model) . LaMEM is a 3-D staggered-grid finite-difference code, which is based upon PETSc (Portable, Extensible Toolkit for Scientific Computation) . The software is massively parallel and thus optimized for the use of high-performance computing (HPC) clusters to enable the computation of high-resolution models in a reasonable amount of time. Here, the matrix is considered impermeable and constrained to no-flow conditions by forcing all three velocity components to be zero at the matrix–void boundary. Besides, the staggered-grid discretization scheme is rescaled at the fluid–matrix interface to provide higher accuracy . Different pressures are applied on two opposing boundaries ($\mathrm{\Delta }P=\mathrm{0.01}$ Pa for all models), while the remaining boundaries are set to no slip. This fixes the principal direction of fluid movement (here, it is in the $y$ direction, e.g., Fig. ). After ensuring that the numerically converged solution is obtained see Appendix A in, the velocity component parallel to the principal flow direction is integrated over the volume to compute the volume average velocity $\stackrel{\mathrm{‾}}{v}$ according to

12$$\stackrel{\mathrm{‾}}{v}={\displaystyle \frac{\mathrm{1}}{V}}\underset{V}{\int }\left|v_y\right|\mathrm{d}y,$$ with the domain volume $V$. To finally obtain the intrinsic permeability $k_i$, $\stackrel{\mathrm{‾}}{v}$ is substituted into Darcy's law for flow through porous media, similar to the approach of : 13$$k_i={\displaystyle \frac{\mathit{\mu }\stackrel{\mathrm{‾}}{v}}{\mathrm{\Delta }P}},$$ with the fluid's dynamic viscosity $\mathit{\mu }$.

Figure 3

Model setup employed in the numerical simulations. The two boundaries where the pressure gradient is applied are indicated. The green velocity vectors are used for the computation of $\stackrel{\mathrm{‾}}{v}$ and scaled according to their magnitude. The subfigure illustrates the staggered-grid discretization scheme of a single voxel.

[Figure omitted. See PDF]

2.3 Synthetic fracture dataset

As in situ data of fractures are rarely accessible and limited to the size of drill cores, numerical studies commonly rely on a stochastic generation approach for synthetic fractures e.g.,. Here, we numerically generate isotropic self-affine surfaces with a MATLAB script , which makes use of random theory and fractal modeling techniques . It uses the standard deviation of surface heights ${\mathit{\sigma }}_{\mathrm{s}}$, the Hurst exponent $H$, the cutoff length $l_{\mathrm{c}}$, and the system size $L$ in $x$ and $y$ directions as input parameters to obtain $s_{\mathrm{u}}(x,y)$ and $s_{\mathrm{l}}(x,y)$, which are built from two independent Gaussian random number fields. The code is slightly modified, such that the seeds for the random number generator are prescribed to produce reproducible results. The mean planes ($x$–$y$ coordinate plane in both cases) of $s_{\mathrm{u}}$ and $s_{\mathrm{l}}$ are separated by varying values of $d$ according to Eq. () to simulate different closure stages of the fracture. Since ignoring mechanical deformation is a common practice , we also assume that both surfaces are in contact at the locations where they overlap. Finally, the data are transferred into a 3-D voxel space of $\mathrm{512}\times \mathrm{512}\times \mathrm{128}$ voxels with a fixed physical voxel size of $\mathrm{0.1}$ mm, resulting in a model domain of $\mathrm{51.2}\times \mathrm{51.2}\times \mathrm{12.8}$ mm. The relative closure and the effective surface area are computed according to Eqs. () and (), respectively. It is important to note that both quantities are computed only within the effective pore space of the fracture parallel to the direction of the applied pressure gradient because only this contributes to the overall flow. For some configurations, it might be possible to have small amounts of trapped pore space within the fracture, which must be excluded before further numerical treatment. We separate the datasets generated in this study into two groups, each with specific sets of input-parameter combinations and a certain number of realizations (input values are listed in Table ). The first group is used to analyze the percolation probabilities and determine the parameter boundaries for geometries of group 2, which are later used for numerical flow simulations. To check whether each fracture configuration in group 1 is able to transmit fluids, i.e., if it is percolating or not, we apply a recursive flood-filling MATLAB routine e.g.,. Then, the percolation probability $p$ represents the mean value of $n$ fracture realizations built from one specific input parameter combination, such that

14$$p={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n{p}_i\text{, with}{p}_i=\left\{\begin{array}{l}\mathrm{1}\text{if percolating}\\ \mathrm{0}\text{if non-percolating}.\end{array}\right.$$ Figure  shows the percolation probability as a function of relative closure $R$, as there was no notable variation with respect to their effective surface areas. Generally, the percolation probability starts reducing from $\mathrm{1}$ at $R\approx \mathrm{1}$ and converges to zero at $R\approx \mathrm{5}$. Higher $l_{\mathrm{c}}/L$ ratios show earlier convergences to their percolation limit, as visible from the two fitted lines for $l_{\mathrm{c}}/L$ $\mathrm{1}$ and $\mathrm{1}/\mathrm{16}$. From the inset plot, it is evident that the contact fraction of all models only depends on the relative closure $R$; first contact between both walls occurs at $R\ge \mathrm{3}\sqrt{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{a}}$, which is in good accordance with . Following this, we have chosen to limit the fracture geometries for the fluid flow simulations to configurations with $R\le \mathrm{1}$ to (i) exclude non-percolation systems and (ii) limit the effect of the above-mentioned “melting” hypothesis, which intensifies with increasing $R$. To ensure applicability to nature, the input values for group 2 (see Table ) are chosen, such that the resulting fracture geometries are classified from “closed” to “open” joints according to . The resulting parameter ranges for fractures in group 2 can be found in Table . For the numerical fluid flow simulations, we implemented the following workflow. First, we apply a flood-filling algorithm on the initial 3-D model along the $x$ direction. $R$ and $S$ are calculated on the resulting effective pore space, followed by a numerical permeability estimation, as explained above, in the same direction. Then, we rotate the initial model by ${\mathrm{90}}^{\circ }$ in the $x$–$y$ plane, and the procedure explained above is repeated. In this manner, two-directional permeability values for every fracture are obtained, resulting in a total sum of 12 800 fluid-flow simulations.

Table 1

Minimal and maximal input values for parameters $d$, ${\mathit{\sigma }}_{\mathrm{s}}$, $H$, and $l_{\mathrm{c}}/L$. $n$ denotes the total number of increments, including minimum and maximum. Subscripts g1 and g2 indicate data for group 1 or 2, respectively. Thus, multiplying the $n$ values of each parameter gives the total number of parameter combinations ($\mathrm{6000}$ for group 1, $\mathrm{320}$ for group 2). The number of realizations for a set (i.e., the number of different random number seeds used to generate the surfaces with one peculiar parameter combination) are given in the footnotes, resulting in a total of 600 000 and 6400 fracture configurations for groups $\mathrm{1}$ and $\mathrm{2}$, respectively.

Parameter	Dimension	Min${}_{\mathrm{g}\mathrm{1}}$	Max${}_{\mathrm{g}\mathrm{1}}$	$n_{\mathrm{g}\mathrm{1}}$	Min${}_{\mathrm{g}\mathrm{2}}$	Max${}_{\mathrm{g}\mathrm{2}}$	$n_{\mathrm{g}\mathrm{2}}$
$d$	mm	0.01	1	20	0.2	5	4
${\mathit{\sigma }}_{\mathrm{s}}$	mm	0.1	0.6	6	0.2	0.5	4
$H$	–	0.1	1.0	10	0.1	1.0	4
$l_{\mathrm{c}}/L$	–	$\mathrm{1}/\mathrm{16}$	1	5	$\mathrm{1}/\mathrm{16}$	1	5

$r_{\mathrm{g}\mathrm{1}}=\mathrm{100}$; $r_{\mathrm{g}\mathrm{2}}=\mathrm{20}$.

Table 2

Resulting minimal and maximal values for mean aperture ($\stackrel{\mathrm{‾}}{a}$), standard deviation of the aperture field (${\mathit{\sigma }}_{\mathrm{a}}$), contact fraction ($c$), relative fracture closure ($R$), effective surface area ($S$), and numerical fracture permeability ($k_{\mathrm{m}}$) for the fracture geometries in group 2.

Parameter	$\stackrel{\mathrm{‾}}{a}$	${\mathit{\sigma }}_{\mathrm{a}}$	$c$	$R$	$S$	$k_{\mathrm{m}}$
Dimension	m	m	–	–	–	m${}^{\mathrm{2}}$
Min	$\mathrm{1.91}\times {\mathrm{10}}^{-\mathrm{4}}$	$\mathrm{1.16}\times {\mathrm{10}}^{-\mathrm{4}}$	$\mathrm{0}$	$\mathrm{0.03}$	$\mathrm{1.04}$	$\mathrm{6.78}\times {\mathrm{10}}^{-\mathrm{14}}$
Max	$\mathrm{4.96}\times {\mathrm{10}}^{-\mathrm{3}}$	$\mathrm{8.51}\times {\mathrm{10}}^{-\mathrm{4}}$	$\mathrm{0.44}$	$\mathrm{0.99}$	$\mathrm{2.49}$	$\mathrm{7.94}\times {\mathrm{10}}^{-\mathrm{7}}$

Figure 4

Percolation probability $p$ (Eq. ) and the mean fractional amount of contact $c$ as a function of relative closure $R$ for all fracture realization sets of groups 1 and 2. Different shades of blue indicate different $l_{\mathrm{c}}/L$ ratios as given in the legend, whereas black lines show best fits to the data.

[Figure omitted. See PDF]

3 Results

3.1 Hydraulic efficiency

In the following section, we present the results of the numerical fluid flow experiments within the geometries of group 2. For this, we normalize the numerically computed permeabilities ($k_{\mathrm{m}}$) by the permeability predicted by the cubic law ($k_{\mathrm{cl}}$) with equivalent mean aperture $\stackrel{\mathrm{‾}}{a}$ of the associated effective pore space:

15$$\mathit{\chi }={\displaystyle \frac{k_{\mathrm{m}}}{k_{\mathrm{cl}}}}.$$

Consequently, one can use the hydraulic efficiency $\mathit{\chi }$ as the correction factor in Eq. () to apply the cubic law to rough fractures. In that way, a fracture whose configuration is close to the parallel plates geometry shows excellent hydraulic efficiency with $\mathit{\chi }$ close to 1. In the $R$–$S$ space, the parallel plate fracture configuration exclusively corresponds to a single point with coordinates $(\mathrm{0},\mathrm{1})$. Perfect hydraulic efficiencies ($\mathit{\chi }=\mathrm{1}$) were validated by flow simulations in parallel-plate fractures.

The key result of this study, a model that corrects the cubic law in terms of quantified fracture roughness, is proposed in Fig. . Due to the complexity of the results, we fitted a regularized surface with a MATLAB function called “gridfit” to approximate the solution in $R$–$S$ space. The function can interpolate scattered data within a prescribed bounding box and a certain amount of smoothing, resulting in a clearer solution image compared to conventional interpolation techniques. Furthermore, to enable an adequate basis for the fitting, the data were cropped above $R=\mathrm{1.0}$ and $S=\mathrm{2.5}$ to provide sufficient data density in $R$–$S$ space, which reduces the total amount of simulations used for the fitting from 12 800 to 10 292.

Figure 5

The distribution of the hydraulic efficiency $\mathit{\chi }$ for different $l_{\mathrm{c}}/L$ ratios as a function of $R$ and $S$. Both axes limits in panel (a) correspond to panels (b)–(e). Dark blue color indicates poor hydraulic efficiency, whereas lighter color shows increasing accordance with the cubic law. The black contour lines indicate the absolute residuals to the fitted surface (compare with Fig. ).

[Figure omitted. See PDF]

The results display significant deviations from the cubic law approximation, even for fractures where both surfaces are not in contact (i.e., $R\le \mathrm{0.23}$). We obtain the lowest $\mathit{\chi }$ values in regions of high $R$ and $S$, with a general trend of increasing $\mathit{\chi }$ for larger $l_{\mathrm{c}}/L$ ratios. To investigate the hydraulic efficiency fluctuations for similar fracture configurations, the absolute residuals of the fitted surface from Fig.  to the simulated data are computed. As before, we fit a regularized surface through the scattered points from which we extract the displayed contour lines (Fig. ). They indicate that the goodness of the fit reduces with increasing $l_{\mathrm{c}}/L$ ratios, especially in the lower right corner. Hence, the non-uniqueness of the data reduces for lower $l_{\mathrm{c}}/L$ ratios, which can be even better seen in Fig. . A regional maximum residual of about $\mathrm{0.2}$ for fractures with a $l_{\mathrm{c}}/L$ ratio of $\mathrm{1}/\mathrm{16}$ enables a more or less unique parametrization refinement, which is not the case for higher $l_{\mathrm{c}}/L$ ratios. Our finding that low $l_{\mathrm{c}}/L$ ratios show lower reduced variability is consistent with the results of .

Figure 6

Absolute residual values of the fitted surfaces in Fig.  for different $l_{\mathrm{c}}/L$ ratios as a function of $R$ and $S$, binned into equally sized boxes. All axes limits in panel (a) correspond to the ones in panels (b)–(e). Orange boxes indicate the mean absolute residual value of the specific bin, whereas the smaller black boxes on top give the maximum absolute residual.

[Figure omitted. See PDF]

An easy integration of this parametrization refinement into a DFN framework requires a mathematical approximation of the fitted surface shown in Fig. a), which was found by the following equation: 16$$\begin{array}{rl}\mathit{\chi }& =\mathrm{1}-\left(\mathrm{0.4809}\tanh \right(\mathrm{0.5139}S)+\mathrm{0.5408})\\ & \cdot \tanh \left({\displaystyle \frac{R}{\mathrm{39.28}\tanh (-\mathrm{2.451}S)+\mathrm{39.47}}}\right).\end{array}$$ To predict single-fracture permeability, it is only necessary to know the mean and standard deviation of the aperture field ($\stackrel{\mathrm{‾}}{a}$ and ${\mathit{\sigma }}_{\mathrm{a}}$), the fractional amount of surface contact ($c$), and the surface area protruding into the void space (${\mathrm{sa}}_{\mathrm{f}}$). From these values, $R$ and $S$ are computed to infer the hydraulic efficiency $\mathit{\chi }$ with Eq. (), which is then multiplied by the permeability predicted by the cubic law with aperture $\stackrel{\mathrm{‾}}{a}$ (see Eq. ). Figure  demonstrates the accuracy of Eq. () to predict hydraulic efficiencies and accompanying permeabilities for fractures with $l_{\mathrm{c}}/L$ ratios of $\mathrm{1}/\mathrm{16}$. To quantify the hydraulic efficiency fluctuations (${\mathit{\sigma }}_{\mathit{\chi }}$) with respect to its correlation length, we provide a model of the form 17$${\mathit{\sigma }}_{\mathit{\chi }}=(a\times e^R+b)(c\times \tanh (S)+d),$$ with corresponding parameter values given by Table .

Figure 7

Parity plot of predicted versus computed hydraulic efficiency $\mathit{\chi }$ for a total of 2554 fractures with an $l_{\mathrm{c}}/L$ ratio of $\mathrm{1}/\mathrm{16}$. The black line indicates the location of perfect parity. Inlet data give the mean and standard deviation of predicted ($k_{\mathrm{p}}$) over computed ($k_{\mathrm{c}}$) permeabilities that correspond to all data points in the plot.

[Figure omitted. See PDF]

Table 3

Coefficients $a$, $b$, $c$, and $d$ of Eq. (), determined by least-square fitting for fractures with varying $l_{\mathrm{c}}/L$ ratios.

$l_{\mathrm{c}}/L$	$a$	$b$	$c$	$d$
$\mathrm{1}/\mathrm{16}$	0.0428	0.1652	$-\mathrm{0.8226}$	0.8822
$\mathrm{1}/\mathrm{8}$	0.6517	$-\mathrm{0.3135}$	$-\mathrm{0.6751}$	0.6625
$\mathrm{1}/\mathrm{4}$	0.9509	$-\mathrm{0.7343}$	$-\mathrm{0.7852}$	0.7672
$\mathrm{1}/\mathrm{2}$	1.0491	$-\mathrm{0.9632}$	$-\mathrm{1.2065}$	1.1752
1	1.3267	$-\mathrm{1.3174}$	$-\mathrm{1.6613}$	1.6178

3.2 Accuracy of the numerical solution

Numerical inaccuracies in solving the Stokes flow equations related to the resolution of the numerical models potentially have an important impact on the results shown here. For numerical permeability estimations of single fractures, the resolution perpendicular to the aperture field is the most crucial part. As the most relevant roughness features are expressed within the uncorrelated region of a fracture (i.e., where $l_{\mathrm{c}}=L$), it is necessary to examine the numerical error introduced due to resolution loss therein. For that, eight fractures with the size of $\mathrm{4096}\times \mathrm{4096}\times \mathrm{512}$ voxels and a $l_{\mathrm{c}}/L$ ratio of $\mathrm{1}/\mathrm{16}$ are generated in the same manner as explained in Sect. 2.3. For each fracture, 16 subsets are drawn that focus uncorrelated regions of the fracture, resulting in subsets of $\mathrm{256}\times \mathrm{256}\times \mathrm{512}$ voxels. By this, the fracture part oriented perpendicular to the applied pressure gradient is over-resolved by a factor of 2. The resolution of these initial models is then consecutively reduced down to $\mathrm{16}\times \mathrm{16}\times \mathrm{32}$ voxels (see inlay of Fig.  for a workflow sketch) while maintaining a constant dimensional aspect ratio. The resulting permeability at every stage ($k_{\mathrm{r}}$) is then compared to the result at maximal resolution ($k_{max}$), assuming that this represents the most accurate solution. Finally, we compute the error norm according to

18$$\left|\right|{\mathit{\delta }}_k\left|\right|=\left|{\displaystyle \frac{k_{\mathrm{r}}-k_{max}}{k_{max}}}\right|.$$ Ideally, the error norm should get negligible at the highest resolution. Figure  shows the mean error norm of a total of 128 uncorrelated fracture subsets as a function of voxel size $\mathrm{\Delta }r$. A mean error of about 0.01 % at maximal resolution indicates optimal convergence to the most accurate solution, which validates the numerical procedure. The voxel size of the numerical models used in this study ($\mathrm{0.1}$ mm) results in an acceptable mean error of 7.2 %, as indicated in Fig. .

Figure 8

Error norm computed by Eq. () as a function numerical voxel size $\mathrm{\Delta }r$. The orange line depicts the mean error of 128 different fracture subsets discretized with decreasing resolutions. The figure inlay sketches the downsampling procedure: the maximal resolution is consecutively decreased by 16 voxels in the horizontal and 32 voxels in the vertical direction (displayed are the maximal, intermediate, and minimal reduction stages). Dashed gray lines indicate the voxel size and associated error norm of the numerical simulations used to provide the refined single-fracture permeability parameterization. Black lines highlight the convergence rate by showing local slopes as indicated by the attached values.

[Figure omitted. See PDF]

4 Discussion

Many studies report that the cubic law increasingly deviates with increasing relative fracture closure , which is usually attributed to the flow channeling around the contact spots within the fracture, introducing an in-plane tortuosity that reduces the permeability. We quantified the deviations from the cubic law due to vertical roughness features (i.e., amplitudes or Hurst exponents; see Fig. ). The results suggest that with increasing fracture surface area protruding into the fluid phase, more drag force accumulates at the fluid–matrix interface, which resists the flow and leads to reduced permeabilities. It is not possible to capture these vertical variations in the flow field with previous 2-D modeling approaches e.g.,, which then results in a biased prediction and the need for more parameters to ensure an adequate quantification. Figure  highlights this issue by computing the norm $\left|\right|{\mathit{\delta }}_k\left|\right|$ between measured and predicted permeabilities for all fractures in this study. With a mean error of $\mathrm{26.7}$ %, Eq. () delivers a better prediction compared to the mentioned studies.

Figure 9

Mean error norm ${\mathit{\delta }}_k$ (see Eq. ) of all fractures considered in this study for different prediction models. The mean error norm recorded for this study is $\mathrm{0.267}$. Inlet plot shows box and whisker plots incorporating all outliers, i.e., representing minimum and maximum recorded values.

[Figure omitted. See PDF]

Figure 10

A sketch of the three different closure regimes, indicated as a function of $R$ with corresponding $p$ and $c$. The lower part of the figure shows three examples of fluid flow simulations from the indicated regimes. The gray shaded area is fracture void space, whereas the white regions indicate contact between both surfaces. Blue lines depict chosen streamlines, approximated from the resulting Stokes velocity vectors.

[Figure omitted. See PDF]

As already shown by , the uncertainty for predicting fracture flow is a function of its $l_{\mathrm{c}}/L$ ratio. Considering flow predictions for uncorrelated fractures (i.e., $l_{\mathrm{c}}/L\ge \mathrm{1}$) is problematic. Blocked pathways connected to the early appearance of the percolation limit (see Fig. ) or flow-enhancing configurations ($\mathit{\chi }>\mathrm{1}$) as also observed by are producing substantial variations in their hydraulic efficiencies. With decreasing $l_{\mathrm{c}}/L$ ratios, the impact of vertical flow tortuosity on its permeability increases relative to the impact of in-plane tortuosity, as both start to act at comparable scales and generally the fractures exhibit larger portions of flow-inhibiting regions compared to flow-enhancing ones see. On the contrary, the fluctuations in the average flow behavior decrease significantly with decreasing $l_{\mathrm{c}}/L$ ratios. This suggests that predicting hydraulic properties is constrained to fractures whose sizes are significantly greater than their correlation lengths. Theoretically, the correlation length is mainly controlled by shear offset and respective gouge generation . With the assumption of a perfectly matched fracture ($l_{\mathrm{c}}=\mathrm{0}$) at its nucleation stage, it is tempting to propose that most natural fractures actually meet the conditions of low $l_{\mathrm{c}}/L$ ratios and subsequently enable the prediction of their hydraulic properties. However, so far, little is known about naturally existing correlation lengths in fractures, as the imaging of in situ fractures is limited to the size of drill cores. Only report measurements of $l_{\mathrm{c}}$ by analyzing the power spectral densities of composite topographies for two matched profiles on opposing joint surfaces, shedding some light onto their natural ranges. From a mechanical perspective, correlation lengths that are equal to the size of the fracture seem rather unrealistic, considering that the shear displacement $d_{\mathrm{s}}$ of fractures scales with their length $L_{\mathrm{f}}$ according to $d_{\mathrm{s}}=\mathit{\alpha }{L_{\mathrm{f}}}^{\mathrm{0.5}}$ . Using $\mathit{\alpha }$ values between $\mathrm{0.01}$ and $\mathrm{0.001}$, which is about the range for Moros joints in , results in a maximal $l_{\mathrm{c}}/L_{\mathrm{f}}$ ratio of $\mathrm{0.01}$. All of that illustrates that further research on fractures correlation lengths is required because the presence of it is omnipresent in most relevant studies .

5 Conclusions

To understand the effects of fracture surface roughness on fluid flow, we performed numerical simulations of high-resolution 3-D Stokes flow within fractures for a large synthetic dataset. By consolidating varying asperity amplitudes and roughness scaling within a new quantity, which accounts for the effective increase of surface roughness compared to its parallel plate equivalent, we were able to provide a new way to characterize fracture roughness. By combining the effective surface area with the relative fracture closure, we established a two-parameter characterization scheme that reads similar to a phase diagram. It is utilized to quantify the hydraulic efficiency of single fractures empirically, i.e., the correction factor applied to the current state-of-the-art fracture permeability parametrization (cubic law). Our findings confirm the results of , and highlight that predicting fracture flow is constrained to scales of at least 16 times larger than the correlation length. The hydraulic efficiency as a function of effective surface area and fracture closure is given by Eq. (), its variability with respect to the correlation length is given by Eq. () and Table , whereas an overall numerical error of $\mathrm{7.2}$ % has to be considered. Ultimately, we used the percolation probability and contact fractions to classify three different closure regimes (see Fig.  for a sketch) that differ in terms of their hydraulic interpretation:

• i.

  The open regime defines fractures whose surface walls are not in contact with each other (e.g., unconfined dilatant or karstified fractures). In this regime, we generally observe a good agreement of the cubic law with hydraulic efficiency between 70 % and 100 %, and only extreme roughness configurations (e.g., needle-shaped mineral coatings) result in larger deviations.

• ii.

  The contact regime is characterized by fractures exhibiting a rapidly decreasing hydraulic efficiency from 70 % to 10 %, up to 1 % in extreme cases, which is caused by strong three-dimensional channeling due to surface roughness and increasing fracture closure. Likely, this regime is most suitable for subsurface conditions, as a certain amount of contact between both fracture surfaces is required to withstand confining pressures.

• iii.

  The percolation regime incorporates fracture configurations that do not percolate at all due to blocked fluid pathways. Here, we do not incorporate fluid flow data, but it is plausible that the hydraulic efficiency is very poor with a maximum of 25 %, which will quickly converge to 0 % due to the effect of decreasing percolation probability with further closure. We observe the no-flow boundary at $R\approx \mathrm{6}$.

Our results generally help to understand the hydraulic response induced by different types of fracture geometries and refine the parametrization of single-fracture permeability given by the cubic law. Moreover, the developed quantification scheme allows monitoring and parametrizing the hydraulic and geometric evolution of fractures during aperture field-shaping processes. This parametrization can easily be incorporated in a DFN modeling framework to investigate the hydraulic responses at reservoir scales, assuming that the minimal correlation length is no longer than $\mathrm{1}/\mathrm{16}$ of the reservoir size. If DFNs of scales close to the correlation length are considered, fluctuations in the average flow behavior are expected. This can modify network-scale flow connectivity and thus requires additional concepts to compute permeabilities e.g.,.

Code availability

The code is available at https://bitbucket.org/bkaus/lamem/src/master/ (); commit: 9c06e4077439b5492d49d03c27d3a1a5f9b65d32.

Author contributions

MOK wrote the initial draft of the manuscript, performed numerical simulations, analyzed the data, and generated the figures. AAP supervised and designed the study, provided the computational framework, and edited the manuscript. TSB assisted in data fitting and edited the manuscript. BJPK helped design the study, assisted in finding an analytical model, and edited the manuscript.

Competing interests

The authors declare that they have no conflict of interest.

Special issue statement

This article is part of the special issue “Faults, fractures, and fluid flow in the shallow crust”. It is not associated with a conference.

Acknowledgements

This work has been funded by the Federal Ministry of Education and Research (BMBF) program GEO:N, grant no. 03G0865A. The authors gratefully acknowledge the computing time granted on the Mogon II supercomputer at Johannes Gutenberg University Mainz (https://hpc.uni-mainz.de, last access: 27 May 2020). The authors sincerely thank Guido Blöcher and an anonymous referee for reviewing this manuscript.

Financial support

This research has been supported by the Bundesministerium für Bildung und Forschung (grant no. 03G0865A). This open-access publication was funded by Johannes Gutenberg University Mainz.

Review statement

This paper was edited by Randolph Williams and reviewed by Guido Blöcher and one anonymous referee.