Articles | Volume 11, issue 3
Research article
29 May 2020
Research article |  | 29 May 2020

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

Maximilian O. Kottwitz, Anton A. Popov, Tobias S. Baumann, and Boris J. P. Kaus

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 cm2 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, CO2 sequestration, nuclear waste disposal, and groundwater management. Fluid flow in fractured reservoirs is commonly modeled by the discrete-fracture-network (DFN) approach (Bogdanov et al.2003; Klimczak et al.2010; Leung et al.2012; de Dreuzy et al.2012), which relies on knowing the permeability of single fractures. The permeability of a single fracture is often approximated by the well-known cubic law (Snow1969; Witherspoon et al.1980), 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 (Durham et al.2001), erosion (Pyrak-Nolte and Nolte2016), and mineral growth (Kling et al.2017), as well as the surrounding stress field (Zimmerman and Main2004; Azizmohammadi and Matthäi2017) 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 (Patir and Cheng1978; Brown1987; Zimmerman and Bodvarsson1996; Oron and Berkowitz1998) 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 xy 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 (Jin et al.2017; Foroughi et al.2018) or Navier–Stokes (Mourzenko et al.1995; Brush and Thomson2003) simulations revealing the non-uniqueness of previous functional approximations of fracture permeability. Factors such as shear displacement (Kluge et al.2017), tortuosity, and the degree of mismatch between the opposing fracture surfaces (Mourzenko et al.2018) 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 (Boffa et al.1998; Ponson et al.2007). Microfractures in Pomeranian shale featured H of 0.3 and 0.5 (Pluymakers et al.2017), 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 (Power and Tullis1991; Schmittbuhl et al.1995; Bouchaud1997; Renard et al.2006; Candela et al.2012), 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 (Plouraboué et al.1995). However, observations of opposing fracture walls (Brown1995) 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 Méheust and Schmittbuhl (2001, 2003), the ratio between system size (L) and the correlation length (lc) defines whether the fracture has an intrinsic permeability or not. Their statistical approach suggested that permeabilities of uncorrelated fractures (i.e., lc/L=1) are strongly fluctuating and anisotropic for the same roughness configurations, revealing the importance of considering low lcL 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 Zimmerman and Bodvarsson1996), 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. 1) and continuity (Eq. 2) equations, which for steady-state flow conditions are given in compact form by


with the fluid's dynamic viscosity μ, pressure P, and velocity vector v=(vx,vy,vz), , ∇⋅, and 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 = - w a 3 Δ P 12 μ ,

with total volumetric flow rate Q, fracture width w, and pressure gradient along the fracture ΔP (see Zimmerman and Bodvarsson1996 for a more detailed derivation). Combining Eq. (3) with Darcy's law for flow through porous media,

(4) Q = - k A Δ P μ ,

with cross-sectional area A, leaves the intrinsic permeability k of an idealized fracture proportional to its aperture by ka2/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 su(x,y) and sl(x,y) are separated by a constant distance d to form a(x,y) according to

(5) a ( x , y ) = s u ( x , y ) + d 2 - s l ( x , y ) - d 2 .

Depending on d, the surfaces may overlap at certain points and form contact areas, such that a0(x,y) is zero at these locations:

(6) a 0 ( x , y ) = a if  a ( x , y ) > 0 0 if  a ( x , y ) 0 .

Assuming that su and sl are self-affine, the standard deviation σ of the aperture field computed at system size l scales as (Brown1995; Schmittbuhl et al.1995):

(7) σ = β l H if  0 l l c β l c H if  l c l L ,

with the maximal system size L and the correlation length lc. Below lc, the fracture is uncorrelated, and it is well correlated above it (see Fig. 1 for a visual explanation). The prefactor β delivers information about the overall amplitude of the surface roughness. H typically denotes the scaling or Hurst exponent with 0<H1, whereas H=1 corresponds to self-similar and H<1 to self-affine scaling (e.g., Mandelbrot1983). 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=1 mm, σs=0.5 mm, H=0.9) with a lcL ratio of 1∕16. (b) Zoomed-in image showing the uncorrelated part of the aperture field. Units of Lx are in millimeters.

Figure 2Two 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) 0.4 and (b) 0.8. The two statistical parameters (a and σ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 a and σ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).


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 σa by the average aperture field a:

(8) R = σ a a ,

with a defined by

(9) a = 1 L x L y x = 0 L x y = 0 L y a 0 ( x , y ) d x d y .

This quantity or its reciprocal is commonly used to infer the amount of contact between the opposing fracture walls (Patir and Cheng1978; Brown1987; Zimmerman and Bodvarsson1996; Méheust and Schmittbuhl2000). Theoretically, it falls in the range 0<R, whereas R=0 shows perfect accordance with parallel plates and the surfaces are in contact if R(32)-1 (see Brown1987).

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. 2 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 saf to twice the area of its projection on the fracture plane (i.e., 2 times the base area perpendicular to the flow direction) sac and normalize it with the fractional amount of the aperture field that has opened, i.e.,

(10) S = sa f sa c 1 1 - c ,

with c being the contact fraction of the aperture field (i.e., where a0(x,y)=0), which leaves 1S, with S=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., Witherspoon et al.1980) by applying a correction factor according to

(11) k χ a 2 12 .

The approximation of χ in terms of quantified fracture roughness, i.e., χ(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. 1) and continuity (Eq. 2) 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) (Kaus et al.2016). LaMEM is a 3-D staggered-grid finite-difference code, which is based upon PETSc (Portable, Extensible Toolkit for Scientific Computation) (Balay et al.2018). 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 (Eichheimer et al.2019). Different pressures are applied on two opposing boundaries (ΔP=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. 3). After ensuring that the numerically converged solution is obtained (see Appendix A in Eichheimer et al.2019), the velocity component parallel to the principal flow direction is integrated over the volume to compute the volume average velocity v according to

(12) v = 1 V V | v y | d y ,

with the domain volume V. To finally obtain the intrinsic permeability ki, v is substituted into Darcy's law for flow through porous media, similar to the approach of Osorno et al. (2015):

(13) k i = μ v Δ P ,

with the fluid's dynamic viscosity μ.

Figure 3Model 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 v and scaled according to their magnitude. The subfigure illustrates the staggered-grid discretization scheme of a single voxel.


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., Brown1995; Candela et al.2009). Here, we numerically generate isotropic self-affine surfaces with a MATLAB script (Kanafi2016), which makes use of random theory and fractal modeling techniques (Persson et al.2004). It uses the standard deviation of surface heights σs, the Hurst exponent H, the cutoff length lc, and the system size L in x and y directions as input parameters to obtain su(x,y) and sl(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 (xy coordinate plane in both cases) of su and sl are separated by varying values of d according to Eq. (5) to simulate different closure stages of the fracture. Since ignoring mechanical deformation is a common practice (Brown1987; Méheust and Schmittbuhl2000, 2003; Mourzenko et al.2018), 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 512×512×128 voxels with a fixed physical voxel size of 0.1 mm, resulting in a model domain of 51.2×51.2×12.8 mm. The relative closure and the effective surface area are computed according to Eqs. (8) and (10), 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 1). 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., Torbert2016). Then, the percolation probability p represents the mean value of n fracture realizations built from one specific input parameter combination, such that

(14) p = 1 n i = 1 n p i , with p i = 1 if percolating 0 if non-percolating .

Figure 4 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 1 at R≈1 and converges to zero at R≈5. Higher lcL ratios show earlier convergences to their percolation limit, as visible from the two fitted lines for lcL 1 and 1∕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 R32σa, which is in good accordance with Brown (1987). Following this, we have chosen to limit the fracture geometries for the fluid flow simulations to configurations with R≤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 1) are chosen, such that the resulting fracture geometries are classified from “closed” to “open” joints according to Bieniawski (1989). The resulting parameter ranges for fractures in group 2 can be found in Table 2. 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 90 in the xy 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 1Minimal and maximal input values for parameters d, σs, H, and lcL. 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 (6000 for group 1, 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 1 and 2, respectively.

rg1=100; rg2=20.

Download Print Version | Download XLSX

Table 2Resulting minimal and maximal values for mean aperture (a), standard deviation of the aperture field (σa), contact fraction (c), relative fracture closure (R), effective surface area (S), and numerical fracture permeability (km) for the fracture geometries in group 2.

Download Print Version | Download XLSX

Figure 4Percolation probability p (Eq. 14) 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 lcL ratios as given in the legend, whereas black lines show best fits to the data.


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 (km) by the permeability predicted by the cubic law (kcl) with equivalent mean aperture a of the associated effective pore space:

(15) χ = k m k cl .

Consequently, one can use the hydraulic efficiency χ as the correction factor in Eq. (11) 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 χ close to 1. In the RS space, the parallel plate fracture configuration exclusively corresponds to a single point with coordinates (0,1). Perfect hydraulic efficiencies (χ=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. 5. Due to the complexity of the results, we fitted a regularized surface with a MATLAB function called “gridfit” (D'Errico2006) to approximate the solution in RS 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=1.0 and S=2.5 to provide sufficient data density in RS space, which reduces the total amount of simulations used for the fitting from 12 800 to 10 292.

Figure 5The distribution of the hydraulic efficiency χ for different lcL 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. 6).


The results display significant deviations from the cubic law approximation, even for fractures where both surfaces are not in contact (i.e., R≤0.23). We obtain the lowest χ values in regions of high R and S, with a general trend of increasing χ for larger lcL ratios. To investigate the hydraulic efficiency fluctuations for similar fracture configurations, the absolute residuals of the fitted surface from Fig. 5 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. 5). They indicate that the goodness of the fit reduces with increasing lcL ratios, especially in the lower right corner. Hence, the non-uniqueness of the data reduces for lower lcL ratios, which can be even better seen in Fig. 6. A regional maximum residual of about 0.2 for fractures with a lcL ratio of 1∕16 enables a more or less unique parametrization refinement, which is not the case for higher lcL ratios. Our finding that low lcL ratios show lower reduced variability is consistent with the results of Méheust and Schmittbuhl (2003).

Figure 6Absolute residual values of the fitted surfaces in Fig. 5 for different lcL 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.


An easy integration of this parametrization refinement into a DFN framework requires a mathematical approximation of the fitted surface shown in Fig. 5a), which was found by the following equation:

(16) χ = 1 - ( 0.4809 tanh ( 0.5139 S ) + 0.5408 ) tanh R 39.28 tanh ( - 2.451 S ) + 39.47 .

To predict single-fracture permeability, it is only necessary to know the mean and standard deviation of the aperture field (a and σa), the fractional amount of surface contact (c), and the surface area protruding into the void space (saf). From these values, R and S are computed to infer the hydraulic efficiency χ with Eq. (16), which is then multiplied by the permeability predicted by the cubic law with aperture a (see Eq. 11). Figure 7 demonstrates the accuracy of Eq. (16) to predict hydraulic efficiencies and accompanying permeabilities for fractures with lcL ratios of 1∕16. To quantify the hydraulic efficiency fluctuations (σχ) with respect to its correlation length, we provide a model of the form

(17) σ χ = ( a × e R + b ) ( c × tanh ( S ) + d ) ,

with corresponding parameter values given by Table 3.

Figure 7Parity plot of predicted versus computed hydraulic efficiency χ for a total of 2554 fractures with an lcL ratio of 1∕16. The black line indicates the location of perfect parity. Inlet data give the mean and standard deviation of predicted (kp) over computed (kc) permeabilities that correspond to all data points in the plot.


Table 3Coefficients a, b, c, and d of Eq. (17), determined by least-square fitting for fractures with varying lcL ratios.

Download Print Version | Download XLSX

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 lc=L), it is necessary to examine the numerical error introduced due to resolution loss therein. For that, eight fractures with the size of 4096×4096×512 voxels and a lcL ratio of 1∕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 256×256×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 16×16×32 voxels (see inlay of Fig. 8 for a workflow sketch) while maintaining a constant dimensional aspect ratio. The resulting permeability at every stage (kr) is then compared to the result at maximal resolution (kmax), assuming that this represents the most accurate solution. Finally, we compute the error norm according to

(18) | | δ k | | = k r - k max k max .

Ideally, the error norm should get negligible at the highest resolution. Figure 8 shows the mean error norm of a total of 128 uncorrelated fracture subsets as a function of voxel size Δ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 (0.1 mm) results in an acceptable mean error of 7.2 %, as indicated in Fig. 8.

Figure 8Error norm computed by Eq. (18) as a function numerical voxel size Δ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.


4 Discussion

Many studies report that the cubic law increasingly deviates with increasing relative fracture closure (Patir and Cheng1978; Brown1987; Zimmerman and Bodvarsson1996), 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. 5). 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., Patir and Cheng1978; Brown1987; Renshaw1995; Zimmerman and Bodvarsson1996), which then results in a biased prediction and the need for more parameters to ensure an adequate quantification. Figure 9 highlights this issue by computing the norm ||δk|| between measured and predicted permeabilities for all fractures in this study. With a mean error of 26.7 %, Eq. (16) delivers a better prediction compared to the mentioned studies.

Figure 9Mean error norm δk (see Eq. 18) of all fractures considered in this study for different prediction models. The mean error norm recorded for this study is 0.267. Inlet plot shows box and whisker plots incorporating all outliers, i.e., representing minimum and maximum recorded values.


Figure 10A 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.


As already shown by Méheust and Schmittbuhl (2003), the uncertainty for predicting fracture flow is a function of its lcL ratio. Considering flow predictions for uncorrelated fractures (i.e., lc/L1) is problematic. Blocked pathways connected to the early appearance of the percolation limit (see Fig. 4) or flow-enhancing configurations (χ>1) as also observed by Méheust and Schmittbuhl (2000) are producing substantial variations in their hydraulic efficiencies. With decreasing lcL 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 Méheust and Schmittbuhl2000). On the contrary, the fluctuations in the average flow behavior decrease significantly with decreasing lcL 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 (Brown1995; Méheust and Schmittbuhl2000). With the assumption of a perfectly matched fracture (lc=0) at its nucleation stage, it is tempting to propose that most natural fractures actually meet the conditions of low lcL 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 Brown (1995) report measurements of lc 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 ds of fractures scales with their length Lf according to ds=αLf0.5 (Schultz et al.2008). Using α values between 0.01 and 0.001, which is about the range for Moros joints in Schultz et al. (2008), results in a maximal lcLf ratio of 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 (Brown1995; Mourzenko et al.1996; Méheust and Schmittbuhl2003; de Dreuzy et al.2012; Pyrak-Nolte and Nolte2016; Mourzenko et al.2018).

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 Méheust and Schmittbuhl (2003), 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. (16), its variability with respect to the correlation length is given by Eq. (17) and Table 3, whereas an overall numerical error of 7.2 % has to be considered. Ultimately, we used the percolation probability and contact fractions to classify three different closure regimes (see Fig. 10 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≈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 1∕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., de Dreuzy et al.2012).

Code availability

The code is available at (Kaus2018); 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.


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 (, 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.


Azizmohammadi, S. and Matthäi, S. K.: Is the permeability of naturally fractured rocks scale dependent?, Water Resour. Res., 53, 8041–8063, 2017. a

Balay, S., Abhyankar, S., Adams, M., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Dener, A., Eijkhout, V., Gropp, W. D., Karpeyev, D., Kaushik, D., Knepley, M. G., May, D. A., Curfman McInnes, L., Tran Mills, R., Munson, T., Rupp, K., Sanan, P., Smith, B. F., Zampini, S., Zhang, H., and Zhang, H.: PETSc Users Manual: Revision 3.10, Tech. rep., Argonne National Lab. (ANL), Argonne, IL, USA, 2018. a

Bieniawski, Z. T.: Engineering rock mass classifications: a complete manual for engineers and geologists in mining, civil, and petroleum engineering, John Wiley & Sons, New York, 1989. a

Boffa, J. M., Allain, C., and Hulin, J. P.: Experimental analysis of fracture rugosity in granular and compact rocks, Eur. Phys. J.-Appl. Phys., 2, 281–289, 1998. a

Bogdanov, I. I., Mourzenko, V. V., Thovert, J.-F., and Adler, P. M.: Effective permeability of fractured porous media in steady state flow, Water Resour. Res., 39, 1023,, 2003. a

Bouchaud, E.: Scaling properties of cracks, Journal of Physics: Condensed Matter, 9, 4319–4344, 1997. a

Brown, S. R.: Fluid flow through rock joints: the effect of surface roughness, J. Geophys. Res.-Sol. Ea., 92, 1337–1347, 1987. a, b, c, d, e, f, g

Brown, S. R.: Simple mathematical model of a rough fracture, J. Geophys. Res.-Sol. Ea., 100, 5941–5952, 1995. a, b, c, d, e, f

Brush, D. J. and Thomson, N. R.: Fluid flow in synthetic rough-walled fractures: Navier-Stokes, Stokes, and local cubic law simulations, Water Resour. Res., 39, 1085,, 2003. a

Candela, T., Renard, F., Bouchon, M., Brouste, A., Marsan, D., Schmittbuhl, J., and Voisin, C.: Characterization of fault roughness at various scales: Implications of three-dimensional high resolution topography measurements, Pure Appl. Geophys., 166, 1817–1851, 2009. a

Candela, T., Renard, F., Klinger, Y., Mair, K., Schmittbuhl, J., and Brodsky, E. E.: Roughness of fault surfaces over nine decades of length scales, J. Geophys. Res.-Sol. Ea., 117, B08409,, 2012. a

de Dreuzy, J.-R., Méheust, Y., and Pichot, G.: Influence of fracture scale heterogeneity on the flow properties of three-dimensional discrete fracture networks (DFN), J. Geophys. Res.-Sol. Ea., 117, B11207,, 2012. a, b, c

D'Errico, J.: Surface Fitting using gridfit, MATLAB Central File Exchange, available at: (last access: 27 May 2020), 2006. a

Durham, W. B., Bourcier, W. L., and Burton, E. A.: Direct observation of reactive flow in a single fracture, Water Resour. Res., 37, 1–12, 2001. a

Eichheimer, P., Thielmann, M., Popov, A., Golabek, G. J., Fujita, W., Kottwitz, M. O., and Kaus, B. J. P.: Pore-scale permeability prediction for Newtonian and non-Newtonian fluids, Solid Earth, 10, 1717–1731,, 2019. a, b

Foroughi, S., Jamshidi, S., and Pishvaie, M. R.: New Correlative Models to Improve Prediction of Fracture Permeability and Inertial Resistance Coefficient, Transport Porous Med., 121, 557–584, 2018. a

Jin, Y., Dong, J., Zhang, X., Li, X., and Wu, Y.: Scale and size effects on fluid flow through self-affine rough fractures, Int. J. Heat Mass Tran., 105, 443–451, 2017. a

Kanafi, M. M.: Surface generator: artificial randomly rough surfaces, MATLAB Central File Exchange, available at: (last access: 27 May 2020), 2016. a

Kaus, B.: LaMEM – Lithosphere and Mantle Evolution Model, available at:, last access: 27 April 2018. a

Kaus, B., Popov, A. A., Baumann, T., Pusok, A., Bauville, A., Fernandez, N., and Collignon, M.: Forward and inverse modelling of lithospheric deformation on geological timescales, in: Proceedings of NIC Symposium, Forschungszentrum Jülich GmbH, Jülich, available at: (last access: 27 May 2020), 2016. a

Klimczak, C., Schultz, R. A., Parashar, R., and Reeves, D. M.: Cubic law with aperture-length correlation: implications for network scale fluid flow, Hydrogeol. J., 18, 851–862, 2010. a

Kling, T., Schwarz, J.-O., Wendler, F., Enzmann, F., and Blum, P.: Fracture flow due to hydrothermally induced quartz growth, Adv. Water Resour., 107, 93–107, 2017. a

Kluge, C., Milsch, H., and Blöcher, G.: Permeability of displaced fractures, Enrgy. Proced., 125, 88–97, 2017. a

Leung, C., Hoch, A., and Zimmerman, R.: Comparison of discrete fracture network and equivalent continuum simulations of fluid flow through two-dimensional fracture networks for the DECOVALEX–2011 project, Mineral. Mag., 76, 3179–3190, 2012. a

Mandelbrot, B. B.: The fractal geometry of nature, Vol. 173, WH freeman, New York, 1983. a

Méheust, Y. and Schmittbuhl, J.: Flow enhancement of a rough fracture, Geophys. Res. Lett., 27, 2989–2992, 2000. a, b, c, d, e

Méheust, Y. and Schmittbuhl, J.: Geometrical heterogeneities and permeability anisotropy of rough fractures, J. Geophys. Res.-Sol. Ea., 106, 2089–2102, 2001. a

Méheust, Y. and Schmittbuhl, J.: Scale effects related to flow in rough fractures, Pure Appl. Geophys., 160, 1023–1050, 2003. a, b, c, d, e, f

Mourzenko, V. V., Thovert, J.-F., and Adler, P. M.: Permeability of a single fracture; validity of the Reynolds equation, J. Phys. II, 5, 465–482, 1995. a

Mourzenko, V. V., Thovert, J.-F., and Adler, P. M.: Geometry of simulated fractures, Phys. Rev. E, 53, 5606,, 1996. a

Mourzenko, V. V., Thovert, J.-F., and Adler, P. M.: Conductivity and Transmissivity of a Single Fracture, Transport Porous Med., 123, 235–256, 2018.  a, b, c

Oron, A. P. and Berkowitz, B.: Flow in rock fractures: The local cubic law assumption reexamined, Water Resour. Res., 34, 2811–2825, 1998. a

Osorno, M., Uribe, D., Ruiz, O. E., and Steeb, H.: Finite difference calculations of permeability in large domains in a wide porosity range, Arch. Appl. Mech., 85, 1043–1054, 2015. a

Patir, N. and Cheng, H.: An average flow model for determining effects of three-dimensional roughness on partial hydrodynamic lubrication, J. Lubr. Technol., 100, 12–17, 1978. a, b, c, d

Persson, B., Albohr, O., Tartaglino, U., Volokitin, A., and Tosatti, E.: On the nature of surface roughness with application to contact mechanics, sealing, rubber friction and adhesion, Journal of Physics: Condensed Matter, 17, R1–R62, 2004. a

Plouraboué, F., Kurowski, P., Hulin, J.-P., Roux, S., and Schmittbuhl, J.: Aperture of rough cracks, Phys. Rev. E, 51, 1675–1685, 1995. a

Pluymakers, A., Kobchenko, M., and Renard, F.: How microfracture roughness can be used to distinguish between exhumed cracks and in-situ flow paths in shales, J. Struct. Geol., 94, 87–97, 2017. a

Ponson, L., Auradou, H., Pessel, M., Lazarus, V., and Hulin, J.-P.: Failure mechanisms and surface roughness statistics of fractured Fontainebleau sandstone, Phys. Rev. E, 76, 036108,, 2007. a

Power, W. L. and Tullis, T. E.: Euclidean and fractal models for the description of rock surface roughness, J. Geophys. Res.-Sol. Ea., 96, 415–424, 1991. a

Pyrak-Nolte, L. J. and Nolte, D. D.: Approaching a universal scaling relationship between fracture stiffness and fluid flow, Nat. Commun., 7, 10663,, 2016. a, b

Renard, F., Voisin, C., Marsan, D., and Schmittbuhl, J.: High resolution 3D laser scanner measurements of a strike-slip fault quantify its morphological anisotropy at all scales, Geophys. Res. Lett., 33, L04305,, 2006. a

Renshaw, C. E.: On the relationship between mechanical and hydraulic apertures in rough‐walled fractures. J. Geophys. Res.-Sol. Ea., 100, 24629–24636, 1995. a

Schmittbuhl, J., Schmitt, F., and Scholz, C.: Scaling invariance of crack surfaces, J. Geophys. Res.-Sol. Ea., 100, 5953–5973, 1995. a, b

Schultz, R. A., Soliva, R., Fossen, H., Okubo, C. H., and Reeves, D. M.: Dependence of displacement–length scaling relations for fractures and deformation bands on the volumetric changes across them, J. Struct. Geol., 30, 1405–1411, 2008. a, b

Snow, D. T.: Anisotropie permeability of fractured media, Water Resour. Res., 5, 1273–1289, 1969. a

Torbert, S.: Applied computer science, Springer, Cham, 2016. a

Witherspoon, P. A., Wang, J. S., Iwai, K., and Gale, J. E.: Validity of cubic law for fluid flow in a deformable rock fracture, Water Resour. Res., 16, 1016–1024, 1980. a, b

Zimmerman, R. W. and Bodvarsson, G. S.: Hydraulic conductivity of rock fractures, Transport Porous Med., 23, 1–30, 1996. a, b, c, d, e, f

Zimmerman, R. W. and Main, I.: Hydromechanical behavior of fractured rocks, International Geophysics Series, 89, 363–422, 2004. a

Short summary
In this study, we conducted 3-D numerical simulations of fluid flow in synthetically generated fractures that statistically reflect geometries of naturally occurring fractures. We introduced a non-dimensional characterization scheme to relate fracture permeabilities estimated from the numerical simulations to their geometries in a unique manner. By that, we refined the scaling law for fracture permeability, which can be easily integrated into discrete-fracture-network (DFN) modeling approaches.