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

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 up on statistical variations of the aperture field. However, often they exhibit non-unique solutions, attributed to the geo5 metrical 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 3D fracture voxel models. Each fracture consists of two random, 25cm wide self-affine surfaces with prescribed roughness amplitude, scaling exponent, 10 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 depends on the relative fracture closure. 15 Next, Stokes flow calculations are performed, using a 3D 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 2D solution fields reveal decreasing cubic-law accordance’s 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 behaviour of 20 hydraulic efficiencies and respective fracture permeabilities as a function of their statistical properties. A model to capture fluctuations around that average behaviour 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 25

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.

Introduction
The geometrical inhomogeneities of single fractures and their effect on fluid flow remains a crucial parameter for understanding the hydraulic properties of fractured reservoirs, such as crystalline or tight carbonate rocks with nearly impermeable matrices. 30 Hence, it has wide-ranging industrial applicability in the fields of petroleum and gas production, geothermal energy recovery, CO 2 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 (Snow, 1969;Witherspoon et al., 1980), assuming that a fracture is com-35 posed 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 Nolte, 2016) and mineral growth (Kling et al., 2017) as well as the surrounding stress field (Zimmerman and Main, 2004;Azizmohammadi and Matthäi, 2017) further modify the geometry of a fracture, causing deviations of the parallel plate assumption. 40 Considerable effort has been made to study the effect of fracture surface roughness on flow and reactive transport behavior.
Early attempts (Patir and Cheng, 1978;Brown, 1987;Zimmerman and Bodvarsson, 1996;Oron and Berkowitz, 1998) employed the 2D 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 45 power led to numerical improvements, with 3D Lattice Boltzmann (Jin et al., 2017;Foroughi et al., 2018) or Navier-Stokes (Mourzenko et al., 1995;Brush and Thomson, 2003) 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 quan-50 tifying surface roughness properties from thin-section-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). Micro-fractures 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 Tullis, 1991;Schmittbuhl 55 et al., 1995;Bouchaud, 1997;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 (Brown, 1995) have demonstrated that the two surfaces tend to be well correlated above a specific length scale and noncorrelated below it, which poses an upper limit to the self-affine scaling in nature. Following Méheust and Schmittbuhl (2001, 60 2003), the ratio between system size L, and the correlation length l c defines whether the fracture has an intrinsic permeability or not. Their statistical approach suggested that permeabilities of uncorrelated fractures (i.e., l c /L = 1) are strongly fluctuating and anisotropic for the same roughness configurations, revealing the importance of considering low l 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 65 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 3D 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 parameter-70 ization of single fracture permeability is proposed, which is demonstrated to provide accurate predictions for the permeability of rough fractures.

Fluid flow in self-affine fractures
Generally, the flow of an incompressible Newtonian fluid is most accurately described by the Navier-Stokes equations (NSE).

75
Assuming, that the flow is solely laminar (Reynolds numbers below unity according to Zimmerman and Bodvarsson, 1996), the fluid viscosity is constant and gravity is negligible at the system size, they reduce to the simpler Stokes equations, i.e., momentum balance (1) and continuity (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 = (v x , v y , v z ), ∇, ∇·, and ∇ 2 denote the gradient, divergence, and Laplace operator for 3D 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 85 law: with total volumetric flow rate Q, fracture width w and pressure gradient along the fracture ∆P (see Zimmerman and Bodvarsson (1996) for a more detailed derivation). Combining equation 3 with Darcy's law for flow through porous media: 90 with cross-sectional area A, leaves the intrinsic permeability k of an idealized fracture proportional to its aperture by k ∝ a 2 /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 u (x, y) and s l (x, y) are separated by a constant distance d to form a(x, y) according to: Depending on d, the surfaces may overlap at certain points and form contact areas, such that a 0 (x, y) is zero at these locations: Assuming, that s u and s l are self-affine, the standard deviation σ of the aperture field computed at system size l scales as (Brown, 1995;Schmittbuhl et al., 1995): 100 with the maximal system size L and the correlation length l c . Below l c , 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 < H ≤ 1, whereas H = 1 corresponds to self-similar and H < 1 to self-affine scaling (e.g., Mandelbrot , 1982). 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.

105
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ā: This quantity or its reciprocal is commonly used to infer the amount of contact between the opposing fracture walls (Patir and Cheng, 1978;Brown, 1987;Zimmerman and Bodvarsson, 1996;Méheust and Schmittbuhl, 2000). 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 ≥ (3 √ 2) −1
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 120 area of the fracture sa f to twice the area of its projection on the fracture plane (i.e., two times the base area perpendicular to the flow direction) sa c and normalize it with the fractional amount of the aperture field that has opened, i.e.: with c being the contact fraction of the aperture field (i.e. where a 0 (x, y) = 0), which leaves 1 ≤ S ≤ ∞, with S = 1, showing perfect accordance with parallel plates. To finally quantify the influence of fracture roughness on its intrinsic permeability, the 125 proportionality resulting from the cubic law needs to be corrected (e.g., Witherspoon et al., 1980) by applying a correction factor according to: The approximation of χ in terms of quantified fracture roughness, i.e. χ(R, S), is the main subject of this study.

Numerical permeability estimation 130
To simulate the laminar flow of an incompressible, isothermal, and isoviscous fluid, we use a 3D binary voxel model input. We solve the linearized momentum balance (eq. 1) and continuity (eq. 2) equation in 3D, using velocity and pressure as primary variables. The coupled system is implemented in the open-source software package LaMEM (Kaus et al., 2016). LaMEM is a 3D staggered grid finite-difference code, which is based upon PETSc (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 y-direction, e.g. Fig. 3). After ensuring that the numerically converged 140 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 velocityv according to: with the domain volume V . To finally obtain the intrinsic permeability k i ,v is substituted into Darcy's law for flow through porous media, similar to the approach of Osorno et al. (2015): with the fluid's dynamic viscosity µ.

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., Brown, 1995;Candela et al., 2010). Here, we numerically generate 150 isotropic self-affine surfaces with a MATLAB script (Kanafi, 2016), 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 l c , and the system size L in x and y direction as input parameters to obtain s u (x, y) and s 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 u and s l are 155 separated by varying values of d according to equation 5 to simulate different closure stages of the fracture. Since ignoring mechanical deformation is a common practice (Brown, 1987;Schmittbuhl, 2000, 2003;Mourzenko et al., 2018), we also assume that both surfaces are in contact at the locations where they overlap. Finally, the data is transferred into a 3D 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 equations 8 and 10, respectively. It is im-160 portant 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, that must be excluded before further numerical treatment.
We separate the datasets generated in this study into two groups, each wich specific sets of input-parameter combinations and certain number of realizations (input values are listed in table 1). The first group is used to analyze the percolation probabilities 165 and determining the parameter boundaries for geometries of group 2, which are later used for 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., Torbert, 2016). Then, the percolation probability p represents the mean value of n fracture realizations built from one specific input parameter combination, such that: 170 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 l c /L ratios show earlier convergences to their percolation limit, as visible from the two fitted lines for l c /L 1 Figure 4. Percolation 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 group 1 and 2. Different shades of blue indicate different lc/L ratios as given in the legend, whereas black lines show best fits to the data. 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 R ≥ 3 √ 2σ a , which is in good accordance with Brown (1987). Following this, we have 175 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 180 flood-filling algorithm on the initial 3D model along 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 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 12800 fluid-flow simulations.

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 m ) by the permeability predicted by the cubic law (k cl ) with equivalent mean apertureā of the associated effective pore space: 190 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 one. In the R-S-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.   195 5. Due to the complexity of the results, we fitted a regularized surface with a MATLAB function called "gridfit" (DErrico, 2006) to approximate the solution in R-S space. The function can interpolate scattered data within a prescribed bounding box  Furthermore, to enable an adequate basis for the fitting, the data was cropped above R = 1.0 and S = 2.5 to provide sufficient data density in R-S-space, which reduces the total amount of simulations used for the fitting from 12800 to 10292.

200
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 l c /L 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 205 we extract the displayed contour lines (Fig. 5). They indicate that the goodness of the fit reduces with increasing l c /L ratios, especially in the lower right corner. Hence, the non-uniqueness of the data reduces for lower l c /L ratios, which can be even better seen in Fig. 6. A regional, maximum residual of about 0.2 for fractures with a l c /L ratio of 1/16 enables a more or less unique parametrization refinement, which is not the case for higher l c /L ratios. Our finding that low l c /L ratios show lower reduced variability is consistent with the results of Méheust and Schmittbuhl (2003).

210
An easy integration of this parametrization refinement into a DFN framework requires a mathematical approximation of the Table 3. Resulting minimal and maximal values for mean aperture (ā), 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. fitted surface shown in Fig. 5 a), which was found by the following equation: To predict single fracture permeability, it is only necessary to know the mean and standard deviation of the aperture field (ā and σ a ), the fractional amount of surface contact (c) and the surface area protruding into the void space (sa f ). From these values,

215
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ā (see eq. 11). Fig. 7 demonstrates the accuracy of eq. 16 to predict hydraulic efficiencies and accompanying permeabilities for fractures with l c /L ratios of 1/16. To quantify the hydraulic efficiency fluctuations (σ χ ) with respect to its correlation length, we provide a model of the form: 220 with corresponding parameter values given by table 3.

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 225 the uncorrelated region of a fracture (i.e., where l c = L), it is necessary to examine the numerical error introduced due to resolution loss therein. For that, eight fractures with the size of 4096x4096x512 voxels and a l c /L ratio of 1/16 are generated in the same manner as explained in section 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 two. The resolution of these initial models is then consecutively reduced 230 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 (k 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: Ideally, the error norm should get negligible at the highest resolution. Fig. 8 shows the mean error norm of a total of 128 235 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.

Discussion
Many studies report that the cubic law increasingly deviates with increasing relative fracture closure (Patir and Cheng, 1978;240 Brown, 1987;Zimmerman and Bodvarsson, 1996), 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 2D 245 modeling approaches (e.g., Patir and Cheng, 1978;Brown, 1987;Renshaw, 1995;Zimmerman and Bodvarsson, 1996), which then results in a biased prediction and the need for more parameters to ensure an adequate quantification. Fig. 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.
As already shown by Méheust and Schmittbuhl (2003), the uncertainty for predicting fracture flow is a function of its l c /L 250 ratio. Considering flow predictions for uncorrelated fractures (i.e. l c /L ≥ 1) 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 l 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 255 ones (see Méheust and Schmittbuhl, 2000). On the contrary, the fluctuations in the average flow behaviour decrease significantly with decreasing l 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 (Brown, 1995;Méheust and Schmittbuhl, 2000). With the assumption of a perfectly matched fracture (l c = 0) at its nucleation stage, it is tempting to propose that most natural fractures actually meet the conditions of 260 low l 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 Brown (1995) report measurements of l c by analyzing the power spectral densities of composite topographies for two matched profiles on opposing joint surfaces, shedding some light into 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 s of fractures 265 scales with their length L f according to d s = αL f 0.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 l c /L f 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 (Brown, 1995;Mourzenko et al., 1996;Méheust and Schmittbuhl, 2003;de Dreuzy et al., 2012;Pyrak-Nolte and Nolte, 2016;Mourzenko et al., 2018).

Conclusions
To understand the effects of fracture surface roughness on fluid flow, we performed numerical simulations of high-resolution 3D Stokes-flow within fractures for a large synthetic dataset. By consolidating varying asperity amplitudes and roughness scaling within a new quantity, that accounts for the effective increase of surface roughness compared to its parallel plate equivalent, we 275 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-ofthe-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 280 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 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 285 karstified fractures). In this regime, we generally observe a good agreement of the cubic law with hydraulic efficiency between 70 − 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−10% up to 1 % in extreme cases, which is caused by strong three-dimensional channeling due to surface roughness and increasing 290 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.

295
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 reser-300 voir scales, assuming that the minimal correlation length is no longer than 1/16 of the reservoir size. If DFN's of scales close to the correlation length are considered, fluctuations in the average flow behaviour 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. 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 305 figures. AAP supervised and designed the study, provided the computational framework and edited the manuscript. TSB assisted to data fitting and edited the manuscript. BJPK helped designing the study, assisted to find an analytical model and edited the manuscript.
Competing interests. The authors declare that they have no competing interests