**Research article**
27 Oct 2021

**Research article** | 27 Oct 2021

# Roughness of fracture surfaces in numerical models and laboratory experiments

Steffen Abe and Hagen Deckert

**Steffen Abe and Hagen Deckert**Steffen Abe and Hagen Deckert

- Institute for geothermal resource management, Berlinstr. 107a, 55411 Bingen, Germany

- Institute for geothermal resource management, Berlinstr. 107a, 55411 Bingen, Germany

**Correspondence**: Steffen Abe (s.abe@igem-energie.de)

**Correspondence**: Steffen Abe (s.abe@igem-energie.de)

Received: 27 Apr 2021 – Discussion started: 19 May 2021 – Revised: 01 Sep 2021 – Accepted: 10 Sep 2021 – Published: 27 Oct 2021

We investigate the influence of stress conditions during fracture formation on the geometry and roughness of fracture surfaces. Rough fracture surfaces have been generated in numerical simulations of triaxial deformation experiments using the discrete element method and in a small number of laboratory experiments on limestone and sandstone samples. Digital surface models of the rock samples fractured in the laboratory experiments were produced using high-resolution photogrammetry. The roughness of the surfaces was analyzed in terms of absolute roughness measures such as an estimated joint roughness coefficient (JRC) and in terms of its scaling properties. The results show that all analyzed surfaces are self-affine but with different Hurst exponents between the numerical models and the real rock samples. Results from numerical simulations using a wide range of stress conditions to generate the fracture surfaces show a weak decrease of the Hurst exponents with increasing confining stress and a larger absolute roughness for transversely isotropic stress conditions compared to true triaxial conditions. Other than that, our results suggest that stress conditions have little influence on the surface roughness of newly formed fractures.

- Article
(8533 KB) -
Supplement
(16301 KB) - BibTeX
- EndNote

It is well known that surfaces of faults and fractures in rocks are rough at all scales (Brown and Scholz, 1985; Hobbs, 1993; Power and Durham, 1997; Candela et al., 2012). The roughness of fracture surfaces is important for a range of geological processes such as the mechanical behavior of faults (Okubo and Dietrich, 1984; Griffith et al., 2010; Candela et al., 2011a, b; Angheluta et al., 2011; Ahmadi et al., 2016) or the fluid flow in jointed rock or fault zones (Chen et al., 2000; Watanabe et al., 2008; Bisdom et al., 2016; Briggs et al., 2017; Jin et al., 2017; Zambrano et al., 2019; Kottwitz et al., 2020). However, the processes and parameters controlling the details of the fracture geometry are not fully understood yet.

Roughness can be defined as the deviation of a surface from a plane. The
degree of roughness of a surface can be described in a number of different
ways, ranging from visual, semi-quantitative approaches such as the “joint
roughness coefficient” (JRC) (Barton, 1973; Barton and Choubey, 1977) to fully
quantitative measures derived directly from the geometrical properties of the
surface such as the root mean square of the first deviation (slope) along a
profile *Z*_{2} (Myers, 1962) or the “structure function” (SF) proposed by
Sayles and Thomas (Sayles and Thomas, 1977). It has been shown that those
measures are closely, but not perfectly, correlated to each other
(Tse and Cruden, 1979; Li and Zhang, 2015). A roughness measure of particular
interest due to its possible use in the parametrization of the fluid flow
properties of rock fractures is the “effective surface area S” proposed by
Kottwitz et al. (2020), which can be considered as an extension of the
“areal roughness index” defined by El-Soudani (1978) and therefore a
2-D analog of the “roughness profile index” defined there (*R*_{p}
in Li and Zhang, 2015). A statistical analysis of rough surfaces shows
that they can often be described as self-affine (Turcotte, 1992; Schmittbuhl et al., 1993, 1995; Bouchaud, 1997; Candela et al., 2009, 2012); i.e., they are statistically invariant
under an affine transformation, but not under a global dilation
(Bouchaud, 1997). In that case, the roughness can be described by a
scaling parameter such as a fractal dimension or a Hurst exponent
(Candela et al., 2009) in addition to a geometric roughness measure such as
the root mean square deviation from an average plane at a given scale. While
most of the previously mentioned parameters, i.e., JRC, *Z*_{2}, *R*_{p}
and SF, are measured along profiles across the surface and are therefore
intrinsically direction dependent, the scaling parameters can be calculated
either directionally or direction independent.

Stress boundary conditions are one of the main factors controlling the shape and structure of faults and fractures in brittle rocks (Faulkner et al., 2010). While some experimental studies have investigated the dependence of the roughness of individual fracture surfaces on the stress conditions under which they were generated (Amitrano and Schmittbuhl, 2002), the use of numerical models makes it much easier to systematically study this issue for a wide range of stress parameters, including those which are difficult to access experimentally.

A large number of numerical modeling approaches has been developed to study the evolution and resulting properties of rough cracks, from statistical approaches like fiber bundle models over lattice methods including random fuse networks (RFNs) and random spring networks (RSNs) to standard continuum-based approaches like finite element models (FEMs) (Alava et al., 2006). In this work, we use numerical simulations based on the discrete element method (DEM) (Cundall and Strack, 1979; Donze et al., 1994; Mora and Place, 1994) to systematically study the formation of fracture surfaces under a wide range of stress conditions and to quantify their geometric properties. The focus of the investigation is on the initial geometry of the freshly formed fracture surfaces, i.e., in the case of shear fractures, before significant slip takes place. This means that the results will be mainly applicable to joints and shear fractures with small displacement, both of which are very common structures in brittle rocks. The DEM approach was chosen due to its particular suitability for the numerical simulation of brittle deformation processes (Mair and Abe, 2008; Schöpfer et al., 2009, 2011; Yoon et al., 2012) and the option to run true triaxial deformation experiments where ${\mathit{\sigma}}_{\mathrm{1}}>{\mathit{\sigma}}_{\mathrm{2}}>{\mathit{\sigma}}_{\mathrm{3}}$, which are difficult to perform in the laboratory. In addition, we compare the results from the DEM simulations with data obtained from the photogrammetric analysis of fracture surfaces generated in triaxial compression experiments in the laboratory.

## 2.1 Discrete element method

To simulate the process of rock fracture under an externally applied loading, we are using the discrete element method (Cundall and Strack, 1979; Donze et al., 1994; Mora and Place, 1994). In this approach, a brittle–elastic material is modeled as a collection of spherical particles interacting with their nearest neighbors either by frictional–elastic interactions or by breakable elastic “bonded” interactions. Based on the force-displacement laws implemented in these interactions the forces and moments acting on each particle can be calculated. The resulting translational and rotational accelerations of the particles are then used to calculate particle movements from Newton's equations. For the breakable bonded interactions, a failure criterion is evaluated, and if the failure threshold has been exceeded, the affected bonded interactions are removed and, if the particles involved are still in contact, replaced by a frictional–elastic interactions.

A range of different implementations exist for each of the interaction types,
differing mainly in the details of the force-displacement law and, in the case of
the bonded interactions, the failure criterion. In this work, we are using a
linear force-displacement law for the normal component of the
frictional–elastic interactions and a Coulomb friction law for the tangential
component as described by Cundall and Strack (1979). For the bonded
interaction, we are using the bond model by Wang et al. (2006) which takes
normal, shear, bending and torsional deformation into account. The stiffness
and strength of the bonds are parameterized using the approach of
Weatherley (2011) which calculates normal, shear, bending and
torsional stiffness from the elastic parameters of an assumed bond material,
specifically from Young's modulus *E*_{b} and Poisson's ratio *ν*_{b},
considering cylindrical bonds with a length and diameter controlled by the
radii of the particles they are connecting. A Mohr–Coulomb failure criterion
is used for the bonds based on the strength parameters of the bond material,
i.e., cohesion *C*_{b} and friction angle Φ_{b}.

Because the size of the individual models is important in this work to obtain high-resolution roughness data from the simulated fracture surfaces, we are using the parallel DEM software ESyS-Particle (Abe et al., 2004), which enables the simulation of sufficiently large models.

## 2.2 Surface extraction

The extraction of surface data from the numerical models requires two steps: (1) the identification of the individual fragments of the sample after fracturing (Fig. 1b) and (2) the calculation which groups of the particles contained in each fragment form an individual fracture surface. The fragments of the broken sample are extracted by constructing an undirected graph from the structure of the DEM model such that the particles form the nodes of the graph and the remaining unbroken bonds form the edges in the graph. The fragments can then be extracted by calculating the connected components of that graph (Abe and Mair, 2005). For each fragment larger than ≈10 % of the original model (Fig. 1c), a ray-casting method is used to determine which of the particles are forming the surface of the fragment. In this approach, a set of parallel lines or “rays” with their origin outside the fragment and a specific direction is defined. The first intersection between each line and one of the particles is calculated using the algorithm proposed by Amanatides and Woo (1987). The positions of the calculated intersection points then form the point cloud describing the fragment surface (Fig. 2).

To get a complete coverage of the fragment surface (Fig. 1e),
i.e., to avoid shadowing effects by “overhanging” parts of the fragment
surface, the calculations are performed for multiple view directions of the
rays. The directions from the mass centers of all neighboring fragments to the
mass center of the fragment and directions deviating from those by 30^{∘}
are used. To identify individual fracture surfaces two additional
post-processing steps are performed. The initial outside surfaces of the
intact model are removed by identifying each particle which was part of the
surface of the model in the initial particle packing and removing the
respective intersection points from the point cloud. In the final step, a
calculation is performed for each particle contributing an intersection point
to the surface point cloud to determine which other fragment is closest to
this particle. This information is then used to split the point cloud into
individual chunks, each representing an individual fracture surface
(Fig. 1f). By performing this step for all fragments in the
model, corresponding pairs of surfaces belonging to the same fracture can be
identified.

The 3-D point clouds generated using this method are collections of $(x,y,z)$ coordinates. However, for most further analysis steps, a representation of the surface as height field relative to a plane, i.e., as ${z}^{\prime}({x}^{\prime},{y}^{\prime})$ is needed. To obtain such a representation, a “best-fit” plane for the point cloud is calculated. The location of such a plane is found by calculating the barycenter $\mathit{b}=({x}_{\mathrm{0}},{y}_{\mathrm{0}},{z}_{\mathrm{0}})$ of the point cloud, i.e.,

and its orientation is determined by the two major eigenvectors *e*_{1}
and *e*_{2} of the covariance matrix **C** of the point
cloud. The third eigenvector of the covariance matrix then determines the
normal *n*_{fp}=*e*_{3} of the plane. Using this, the in-plane
coordinates $({x}^{\prime},{y}^{\prime})$ of each point $\mathit{p}=(x,y,z)$ and its perpendicular
distance *z*^{′} from the plane can be calculated as

It should be noted that a surface can only be represented correctly as a height
field in this way if there are no parts of the surface which are
“overhanging” with respect to the normal of the fitted plane, i.e., if there
are no points on the surface with identical $({x}^{\prime},{y}^{\prime})$ but different
*z*^{′}. However, this is generally the case for the fracture surfaces generated
in the numerical models.

## 2.3 Roughness characterization

A roughness measure commonly used in the study of the mechanical behavior of rock surfaces is the JRC defined initially as a parameter relating the shape of a rock joint to its peak shear strength; see Eq. (9) in Barton (1973) or Eq. (2) in Barton and Choubey (1977). Its relation to the geometry of the joint surfaces was only qualitatively defined by assigning JRC values to a set of standard profiles (Barton and Choubey, 1977, Fig. 8). To estimate the JRC of an arbitrary profile from measured geometrical data, a wide range of empirical formulas have been developed in the literature (Li and Zhang, 2015, Table 2). To calculate the approximate JRC of the fracture surfaces generated in the numerical models and the laboratory experiments for three of the 47 equations presented there have been chosen. The subscript of the JRC in the equations below shows the respective number of the equation in Li and Zhang (2015, Table 2).

where *R*_{p} is the “roughness profile index”,
$\mathit{\delta}={R}_{\mathrm{p}}-\mathrm{1}$ the “profile elongation index” and *Z*_{2} the
“root mean square of the first deviation of the profile”, all as defined in
Li and Zhang (2015, Table 1). *R*_{p} is therefore calculated as

and

where *x*_{i} is the abscissa of profile point *i*, *y*_{i} its height above a
mean value, and *N* is the number of sample points. Given that those parameters
are calculated along profiles, the irregular point clouds generated using the
method described in Sect. 2.2 first need to be mapped
to a regular grid.

Self-affine rough surfaces are characterized by the fact that they are statistically invariant under an affine transformation

where *x*,*y* are the “in plane” coordinates of the surface and *z* is the
“height” of the surface above a given mean plane
(Fig. 3a). The exponent *H* is the Hurst exponent or
“roughness index” (Mandelbrot and Van Ness, 1968; Mandelbrot, 1985; Bouchaud, 1997). A range of different method for the calculation of the Hurst
exponent have been described in the literature (Renard et al., 2006; Candela et al., 2009), most of them either based on the evaluation of the power
spectrum of the surface or correlation functions between the heights of points
on the surface depending on their mutual distance. Because the point clouds
describing the surfaces generated by the approach described in
Sect. 2.2 do not form a regular grid, spectral
methods would require an additional interpolation step. Aside from the
additional computational effort required, this might also introduce some
difficult to quantify errors in the calculation of the Hurst exponent
(Kottwitz, 2017).

In this work we therefore use the “height–height correlation function” method as described by Candela et al. (2009). However, in contrast to the description in Candela et al. (2009), the function is not calculated from 1-D profile data but directly from the 2-D surface. The radially averaged height–height correlation function is calculated as the root mean square (rms) averaged height difference of all point pairs within a given distance range, i.e.,

where $\mathrm{\Delta}r=\sqrt{\mathrm{\Delta}{x}^{\mathrm{2}}+\mathrm{\Delta}{y}^{\mathrm{2}}}$ is the “in-plane” distance
between the points in the pair, Δ*z* is the height difference between
the points (Fig. 3b), *w* is the size of the distance bins over
which the height differences are averaged, and *n* is the number of particle
pairs in the respective distance bin. For the calculation of the angular
dependence of the height–height correlation function, the direction between the
two points of the pair is calculated as $\mathit{\varphi}=\mathrm{arctan}\left(\frac{\mathrm{\Delta}y}{\mathrm{\Delta}x}\right)$ and the summation of the height differences is adjusted from 1-D distance
bins in Eq. (11) to 2-D (distance, direction) bins.

where *w*_{r} is the bin size with respect to the in-plane distance of the
points and *w*_{ϕ} is the bin size with respect to the direction from one
point of the pair to the other. Due to the large number of points contained in
the surface point clouds, which is as large *n*_{p}≈500 000
in some cases, and the resulting computational cost if all of the
${n}_{\mathrm{p}}^{\mathrm{2}}/\mathrm{2}$ particle pairs would be taken into account, only a random
sample of 10 000 points (i.e., $\approx \mathrm{5}\times {\mathrm{10}}^{\mathrm{7}}$ point pairs) is used for
each surface. Tests have shown that the reduction in the number of particle
pairs evaluated does not impact the results.

Given that for a self-affine surface, the height–height correlation function
follows a power law, i.e., *c*(Δ*r*)∝Δ*r*^{H}
(Candela et al., 2009), the Hurst exponent *H* can be calculated by fitting
a linear function to the straight part of the log–log plot of *c*(Δ*r*)
vs. Δ*r*. The slope of the linear function is the Hurst exponent. The
Hurst exponent *H* and the fractal dimension *D* of an object are related as
$D=\mathrm{2}-H$ for a 1-D profile (Mandelbrot, 1985) or, more generally, $D=n+\mathrm{1}-H$,
where *n* is the dimension of the object (Yang and Lo, 1997), i.e., *n*=1 for a
profile and *n*=2 for a surface.

A set of numerical simulations was performed to generate fracture surfaces under a wide range of stress conditions (Fig. 4a) and some natural rock samples were fractured in laboratory experiments (Fig. 4b and c). Both numerically and experimentally generated surfaces have been analyzed using the methods described in Sect. 2.3.

## 3.1 Numerical models

To generate a set of model fracture surfaces, a large number of deformation
experiments have been simulated. The set of simulations consists of unconfined
compression (*σ*_{1}>0, ${\mathit{\sigma}}_{\mathrm{2}}={\mathit{\sigma}}_{\mathrm{3}}=\mathrm{0}$), unconfined tension
(*σ*_{3}<0, ${\mathit{\sigma}}_{\mathrm{1}}={\mathit{\sigma}}_{\mathrm{2}}=\mathrm{0}$), standard triaxial compression
(*σ*_{1}>0, ${\mathit{\sigma}}_{\mathrm{2}}={\mathit{\sigma}}_{\mathrm{3}}>\mathrm{0}$) and true triaxial compression
(${\mathit{\sigma}}_{\mathrm{1}}>{\mathit{\sigma}}_{\mathrm{2}}>{\mathit{\sigma}}_{\mathrm{3}}>\mathrm{0}$) experiments. In all compressive models,
*σ*_{1} is parallel to the *y* axis, *σ*_{2} is parallel to the *x* axis,
and *σ*_{3} is parallel to the *z* axis, whereas in the unconfined tensile
models *σ*_{3}, i.e., the extension direction, is parallel to the *y* axis,
*σ*_{1} is parallel to the *x* axis, and *σ*_{2} is parallel to the
*z* axis.

All models are using box-shaped samples with an aspect ratio of $\mathrm{1}:\mathrm{2}:\mathrm{1}$
contained between two servo-controlled plates in the case of the unconfined
compression and tension experiments or six servo-controlled plates for the
standard triaxial and true triaxial experiments. While deformation experiments
in the laboratory usually use cylindrical samples, we decided in favor of
box-shaped samples because they make it much easier to apply the two different
confining stresses in the true triaxial tests. In the tension experiments, the
plates are connected to the boundary particles of the sample by unbreakable
bonds which only induce a force parallel to the normal of the plate but not
perpendicular to it. This means the particle are free to move parallel to the
loading plate, avoiding heterogeneous deformation (“necking”). In the
compressive experiments, both the axial loading plates and, in the confined
models, the plates along the *x* and *z* surfaces of the sample interact with
the boundary particles by frictionless elastic interactions.

In the unconfined experiments (${\mathit{\sigma}}_{\mathrm{2}}={\mathit{\sigma}}_{\mathrm{3}}=\mathrm{0}$), a simple loading
procedure is used, applying a prescribed displacement rate to the plates at
the *y* ends of the model to produce axial shortening or extension. During an
initial phase, the plate speed is ramped up smoothly according to a cosine
function until the chosen speed is reached and then it is held constant for
the main phase of the experiment. In the confined experiments (${\mathit{\sigma}}_{\mathrm{2}}\ge {\mathit{\sigma}}_{\mathrm{3}}>\mathrm{0}$), this loading procedure is preceded by a ramp-up of the stresses
applied to the plates at the *x*- and *z*-sides of the sample until
*σ*_{xx}=*σ*_{2} and *σ*_{zz}=*σ*_{3}. For the smooth ramp-up of
the applied stresses, the same cosine function is used as for the ramp-up of
the axial deformation rate in order to minimize unnecessary vibrations in the
model. During this phase, a stress is also applied to the loading plates at the
*y* ends of the sample such that *σ*_{1}=*σ*_{2}. After a subsequent
“rest” phase where the stress on all plates is held constant for given time
to allow the particle movement introduced by the initial loading to dissipate,
the same axial shortening as in the unconfined compression experiments is
applied. A range of confining stresses from ${\mathit{\sigma}}_{\mathrm{2}}={\mathit{\sigma}}_{\mathrm{3}}=\mathrm{0}$ to
${\mathit{\sigma}}_{\mathrm{2}}={\mathit{\sigma}}_{\mathrm{3}}=\mathrm{15}$ MPa was used for the numerical models in this work. In
order to avoid the effect of abrasion modifying the roughness of the fracture
surfaces after their initial formation, the state of the model immediately
after one or more through-going fractures have formed was used for the
extraction of fracture surfaces described in
Sect. 2.2. That is, there is no, or at least very little,
post-failure slip on those surfaces.

A model size of $\mathrm{55}\times \mathrm{110}\times \mathrm{55}$ model units was chosen with a
particle sizes ranging from *R*_{min}=0.2 to *R*_{max}=1.0,
resulting in ≈950 000 particles for those models
(Fig. 1a). This model size was found in initial tests to
provide a good balance between model resolution and computational cost. For
the construction of the initial particle arrangement for the models, the
insertion-based packing algorithm by Place and Mora (2001) was used. This
algorithm generates dense particle packings having a power-law particle size
distribution with an exponent of approximately −3; i.e., the number of
particles with given radius *r* is roughly proportional to *r*^{−3}.

In all deformation experiments, the final loading plate speed was set to $\approx \mathrm{17}\phantom{\rule{0.125em}{0ex}}\mathrm{cm}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$. This is significantly higher than in real experiments, but using real lab values ($\mathrm{\mu}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}\mathrm{\dots}\mathrm{mm}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{-\mathrm{1}}$) would lead to impractically long computing times because the time step of the calculations is restricted to values of $\mathrm{\Delta}t\le \mathrm{3}\times {\mathrm{10}}^{-\mathrm{8}}$ s due to numerical stability constraints. Tests have shown that the increased velocities do not significantly influence the model results.

The mechanical properties of the DEM material have been calibrated to values
similar to those of a typical sedimentary rock. The target values, Young's
modulus *E*=30 GPa and unconfined compressive strength
UCS=80 MPa, are within the range of sandstone or medium to
high porosity limestone (Zoback, 2007). The failure strength was found to
vary by less than 1 % among samples. These parameters do not provide
a direct match to the mechanical properties of the rocks used in the
laboratory tests (Sect. 3.2), but the important ratio between
failure strength of the material and the confining stress applied in the
laboratory experiments lies well within the range covered by the numerical
models (Fig. 5b). Because the details of the fracture
behaviors of individual samples in DEM models show a well-known dependence on
the initial random particle arrangement (Koyama and Jing, 2007; Abe et al., 2011; Fakhimi and Gharahbagh, 2011), at least five simulations with different realization of the
particle packing have been performed for each parameter set in order to
quantify this variability.

To improve the coverage of the chosen range of stress conditions, data from a related study (Ohagen, 2019) were integrated into the analysis (Fig. 5a). This study was using an identical model setup, except for slightly smaller models with dimensions of $\mathrm{40}\times \mathrm{80}\times \mathrm{40}$ model units (≈360 000 particles) and $\mathrm{50}\times \mathrm{100}\times \mathrm{50}$ model units (≈710 000 particles) compared to the roughly 950 000 particles used in most models in this work.

## 3.2 Laboratory tests

To compare the roughness of fractures created in the DEM models with the roughness of real fractures, we conducted a number of laboratory uniaxial and triaxial deformation experiments. For our study, we used a suite of fine grained, low porosity Upper Jurassic carbonate rock samples and additionally one Lower Triassic sandstone sample, both from Franconia, Germany. Sample size for the experiments were 55×110 mm cylinders. The main goal of the experiments was to produce fractures for given stress conditions which could then be used for roughness analyses.

The sandstone uniaxial compressive strength (UCS) experiment lead to an
typical hourglass fracture pattern, splitting the sample into a small number
of larger fragments, which could be used for further analyses
(Fig. 4c). Unfortunately, for the UCS and most of the triaxial
experiments of the carbonate rocks, the samples disintegrated into a very large
number of very small fragments, leaving no suitable fracture surfaces to
analyze. See Fig. S1 in the Supplement for a typical example. This applied in
particular to the samples loaded with small confining pressures. Only in one
experiment with a confining pressure of 30 MPa post-deformation
fragments were large enough for our planned fracture surface analyses
(Fig. 4b). From the suitable fragments, we constructed a
digital three-dimensional surface model using photogrammetric methods. The
models were built from more than 100 single pictures of the samples from
different perspectives using a 12-megapixel SLR camera and a 40 mm
macro lens. The photos were taken from a distance of 5 to 10 cm
between front lens and the object, which is close to the minimum focus
distance of the lens used. The models were then clipped to the fracture plane
of interest. The remaining surface geometry was exported as 3-D point cloud
data with approximately 2.2 million data points in total, resulting in a point density of
approximately 28 000 points cm^{−2} and an average point distance of 60 µm.

The generated point clouds were then used for roughness analyses of the
fracture surfaces following the approach described in
Sect. 2.3. Besides the creation of fracture surfaces the
deformation experiments were also used to derive typical geomechanical
properties of the carbonate and sandstone samples which were used for
comparison with the DEM models. For the carbonate rocks, a UCS of ≈285 MPa was obtained and ≈85 MPa for the sandstone
sample. For the limestone, a friction coefficient *μ*=0.7 was derived from
experiments with confining pressures ranging between 0 to
30 MPa. Young's modulus was measured at *E*=48 GPa for the
limestone and *E*=12.5 GPa for the sandstone.

## 4.1 Numerical models

Based on the data produced by a total of 131 numerical simulations, the
geometrical properties of 388 fracture surfaces have been analyzed. The
fracture orientations were as expected under the stress conditions. The dip
angle was typically within 25–35^{∘} of *σ*_{1}, i.e., 55–65^{∘}
assuming *σ*_{1} to be vertical. The strike direction of the majority of
the fractures was within ≈10^{∘} of *σ*_{2} in the true
triaxial models (*σ*_{2}≠*σ*_{3}) and more or less randomly
distributed under transverse isotropic stress conditions (*σ*_{2}=*σ*_{3}).

In an initial step, the joint roughness coefficients for a small set of
surfaces were approximated using Eqs. (5)–(7). The
results did show that the resulting JRC values were consistently above the
range defined by Barton and Choubey (1977), i.e., larger than 20, and
therefore also outside the range of validity of the approximation equations in
Li and Zhang (2015). Similarly, the geometric parameters *R*_{p}
(Eq. 8) and *Z*_{2} (Eq. 9) from which the estimated JRC
values were calculated, were outside the applicable ranges given there. While
the roughness produced by the numerical models is therefore outside the range
for which the fitting equations collected by Li and Zhang (2015) were
originally intended, Fig. 1a in their work suggests that Eq. (5)
would be the best option to extend the range of approximate JRCs to the
surface geometries observed here because it provides a particularly good fit
at large values (i.e., *Z*_{2}≈0.35–0.4, JRC≈20). Therefore, Eq. (5) was used to estimate the average JRC for
each of 261 surfaces based on a total of ≈24 400 profiles. The
remaining 127 of the 388 surfaces studied were found to be too small in at
least one of the dimensions to allow the extraction of sufficiently long
profiles. For each surface, profiles were generated in two orthogonal
directions to check for a possible anisotropy of the surface roughness. The
results did show that the mean estimated JRCs for the profiles differs by less
than 10 % between the two direction, which is generally less than the
standard deviation between the profiles within one direction. Plotting the
estimated JRC for the analyzed surfaces against the mean confining stress of
the models (Fig. 6a) shows that there is no clear trend of
JRC vs. confinement, but that models with transversely isotropic confinement
(*σ*_{2}=*σ*_{3}) generally have higher JRC values than models fractured
under true triaxial conditions, i.e., *σ*_{2}≠*σ*_{3}. The directly
calculated geometric roughness measures, i.e., *R*_{p} and *Z*_{2}, show
a very similar pattern (Fig. 6b and c).

The perpendicular distance or “height” of the points of the fracture surfaces above a fitted fit plane is calculated according to Eq. (4). Analysis shows that the heights are normally distributed (Fig. 7) as expected for fracture surfaces (Ponson et al., 2007), allowing a “rms roughness” ${h}_{\text{rms}}=\frac{\mathrm{1}}{n}\sqrt{{\sum}_{n}\left({{z}^{\prime}}^{\mathrm{2}}\right)}$ to be calculated.

Plotting the rms roughness *h*_{rms} of all models against the mean
confining stress (Fig. 8) shows that there is no clear
dependence between the two parameters, except for a difference between
transverse isotropic (*σ*_{2}=*σ*_{3}) and true triaxial (*σ*_{2}≠*σ*_{3}) stress conditions. In the case of the transverse isotropic confinement,
the observed rms roughness ${h}_{\text{rms}}=\mathrm{2.35}\pm \mathrm{0.78}$ model units is
higher than in the case of true triaxial conditions where ${h}_{\text{rms}}=\mathrm{1.63}\pm \mathrm{0.48}$ model units. It can also be observed that the rms roughness of the
models subjected to unconfined extension (blue marker in
Fig. 8) is smaller at ${h}_{\text{rms}}=\mathrm{1.51}\pm \mathrm{0.44}$
model units than that of the models subjected to unconfined compression with
${h}_{\text{rms}}=\mathrm{2.76}\pm \mathrm{0.88}$ model units. This difference, however, is
possibly at least in part an artifact of the different size of the fracture
surfaces between the two model groups. In the tensile case, the fractures tend
to be roughly normal to the extension direction, i.e., the long axis of the
model and their size is therefore restricted to the small cross section of the
model. In contrast, the fracture surfaces in the compressive case tend to be
oriented such that their normal is at an angle of ≈55…60^{∘} to the compression direction and can therefore grow as large as a
plane diagonally across the model, i.e., more than twice the size compared to
the tensile case. Plotting the height–height correlation function
(Eq. 11) of the surfaces in a log–log plot
(Fig. 9) shows a clear linear section which, for most
surfaces analyzed, ranges from Δ*r*_{min}≈1.5–2 model
units, i.e., somewhat more than the maximum particle size, to about half of the
smaller dimension of the surface, which in most cases means Δ*r*_{max}≈20–30 model units. This linear section in the log–log
plot, representing a power-law dependence *c*(Δ*r*)∝Δ*r*^{h}
shows that the surface is indeed self-affine, at least for the range of scales
covered by the linear section.

In order to verify that the observed self-affine structure of the fracture surfaces generated in the numerical model is indeed a result of the fracture process and not an artifact caused by the intrinsic roughness of surfaces in the particle model, a number of “quasi-planar” surfaces were generated in the particle model and their roughness was analyzed. For this purpose, one of the blocks of packed particles used in the DEM simulations of the triaxial tests (Fig. 1a) was cut with an arbitrarily oriented plane; i.e., the particles on one side of the plane were removed. The remaining fragment of the block then underwent the same surface extraction and roughness analysis procedures as the fracture surfaces produced in the deformation experiments. The result (Fig. 10) shows that the height–height correlation function of the cut surface is essentially flat from the particle scale up to the model size. Performing this analysis on multiple cut surfaces did show that this is independent of the orientation of the cut plane and the details of the particle packing. Only the absolute value of the roughness of the cut surfaces depends somewhat on the size range of the particles. Calculating the joint roughness coefficients for the cut surfaces according to Eq. (5) did, as expected, produce non-zero values of the JRC. However, the JRC values for the cut surfaces are in the range of 11.5–12, which is much smaller than the values observed in the fracture surfaces generated in the numerical models (JRC≈23–32, Fig. 6a). It can therefore be assumed that, while there is some contribution of the intrinsic particle-scale roughness to the total roughness of the fracture surface, the self-affine structure of the fracture surfaces as well as the major part of their total roughness is due to the fracture process and not the particle structure of the model as such.

Performing the calculations for all 388 fracture surfaces extracted from the
numerical models produced Hurst exponents ranging from 0.2 to 0.6. To
investigate possible dependencies on the stress conditions under which the
fractures were created, the average Hurst exponents over all surfaces
generated in each set of simulations with identical boundary conditions have
been calculated. The mean value of *H* for the sets varies between 0.3 and
0.45, with the variation between top and bottom quartile within each set
typically in the range of 0.05…0.1. Due to the relatively small number
of data points within each set of models, i.e., between 8 and 28 surfaces, and
the observed asymmetry of the error distribution in some instances, quartiles
have been calculated and plotted (Figs. 11 and
12) instead of standard deviations. Plotting the calculated
Hurst exponents against the mean confining stress $({\mathit{\sigma}}_{\mathrm{2}}+{\mathit{\sigma}}_{\mathrm{3}})/\mathrm{2}$
(Fig. 11) shows a weak trend towards lower Hurst exponents
with increasing confinement. No dependence of the Hurst exponent on the ratio
between the confining stresses ${\mathit{\sigma}}_{\mathrm{3}}/{\mathit{\sigma}}_{\mathrm{2}}$ could be observed
(Fig. 12).

## 4.2 Laboratory tests

To characterize the roughness of the fracture surfaces produced in the
laboratory deformation tests, we examined the photogrammetrically produced
point clouds of the single sample fragments. For each of the limestone and
sandstone samples, one fracture surface was chosen. The maximum sampling area
for the roughness investigation was 14 cm^{2} for the sandstone and
≈25 cm^{2} for the limestone sample. The analyses of the
heights distances of the single points of the point clouds above their fitted
mean planes revealed a normal distribution of the heights. Thus, a calculation
of the rms roughness is justified (Fig. 13). The
height–height correlation functions of these surfaces have a well-defined
linear section in a log–log plot proving a self-affine geometry in a distance
range between ≈0.1 and ≈1 cm, both for the limestone
and sandstone sample (Figs. S2 and S3 in the Supplement). With distances
larger ≈1 cm, a flattening of rms curve can be observed,
marking the upper end of the power-law relationship between the point distance
and the rms height difference. From the linear slope segments of the
correlation functions, similar Hurst exponents could be deduced with *H*=0.66
for the sandstone and *H*=0.69 for the limestone when analyzing the maximum
sampling area on the respective fracture.

To check whether the size of the investigation area on the fracture surfaces
has an effect on calculated Hurst exponents, we analyzed the height–height
correlation functions and Hurst exponents for a suite of different area sizes
(Fig. 14). For the limestone sample, the mean *H* value
results in *H*=0.73 with a standard deviation of 0.08. The sandstone sample
shows a clearly lower mean *H* value of *H*=0.6 and a standard deviation of
0.05. A stronger scatter of Hurst exponents can be observed for the smallest
analyzed sample area size of ≈1 cm^{2}, ranging between
*H*=0.43 and 0.64 for the sandstone surface and between *H*=0.67 and 0.85 with
two outliers of *H*≈0.5 for the limestone surfaces. For these outliers,
a closer investigation of the corresponding rms–distance curves shows that two
different linear sections could be derived, one with a higher Hurst exponent
for smaller distances and one lower Hurst exponent for larger distances.

For both sample surfaces, the JRC was estimated
using the same methods as for the numerical models. The results show that the
estimated JRC is dependent on the sampling resolution, i.e., the number of
sampling points on the profile, specifically that the calculated value of the
JRC is increasing with smaller sampling intervals
(Fig. 15). This is a known effect which is caused by the
dependence of the underlying geometric parameters *R*_{p} and *Z*_{2},
from which the estimated JRC is calculated, on the sampling interval used
(Yu and Vayssade, 1991). It is also to be expected based on the fact that the
analyzed surfaces are self-affine. In that case, the dependence of
*R*_{p} on the sampling interval is directly described by the
“compass dimension” (Mandelbrot, 1985) of the profile. For the fracture
surfaces in the numerical models (Sect. 4.1), a similar
analysis of the resolution dependence of the JRC was not done because of the
lower intrinsic resolution of the point clouds which limits profiles to less
than 100 sample points in most cases.

The empirical equations used for the calculation of JRC from measured
geometric parameters are usually derived based on sampling resolutions between
100 and 400 points per profile (Tse and Cruden, 1979; Yu and Vayssade, 1991; Li and Zhang, 2015). Specifically, the equation used in this work to estimate JRC
from *Z*_{2} (Eq. 5) was derived by Tse and Cruden
(Tse and Cruden, 1979) using a sample interval of 1.27 mm at a profile
length of 25 cm, i.e., slightly less than 200 points. The surfaces
analyzed here have dimensions of about 7 cm×5 cm for
the sandstone and approximately 10 cm×4.5 cm for the
limestone. Therefore, a sampling interval of between 0.25 and 0.5 mm
will produce a similar number of sample points along the profiles. Therefore,
the best estimates for the average JRC of the fracture surfaces produced in
the laboratory experiments are for the sandstone JRC≈9–11 in
the direction parallel to shortening direction in the deformation experiment
and JRC≈11–13 perpendicular to it
(Fig. 15). For the limestone, the estimates are JRC≈6.5–7.5 in the parallel direction and JRC≈16–17 in
the perpendicular direction. In both cases, the JRC shows a clear anisotropy
between the two directions. However, this anisotropy is much larger in the
limestone compared to the sandstone sample.

The results of the analysis of the simulation data (Sect. 4.1) shows that the roughness of the fracture surfaces generated in the numerical models is high compared to natural rock fractures usually considered in the geomechanical literature. In the numerical models, the surfaces show estimated JRC values larger than 23 and in some case exceeding 30, whereas the JRC for natural surfaces was originally only defined for a range up to 20 (Barton, 1973; Barton and Choubey, 1977). In contrast, the natural rock samples analyzed in this work (Sect. 4.2) show JRC values between 6 and 17, which is well within the range defined by Barton (1973).

However, as described in Sect. 4.1, the JRC values for
the numerical model contain a small contribution due to the intrinsic
particle-scale roughness of the model. If we consider that the total roughness of the
surface is the sum of the roughness due to the particle structure of the
surfaces and the roughness due to the actual fracture process, and if we
assume that those contributions are not spatially correlated with each other,
it would be possible to correct the calculated JRC values by removing the
effect of the particle-scale roughness. The parameter *Z*_{2} on which the
calculation of the JRC is based (Eq. 5) is calculated from the rms
of the first derivative of profiles along the surface (Eq. 9). Based
on the assumption that the particle-scale roughness and the fracture-generated
roughness are not spatially correlated, this means that the total *Z*_{2} is the
rms of the *Z*_{2} values of the two parts, and therefore the value *Z*_{2f} of
the fracture-generated roughness can be estimated as
${Z}_{\mathrm{2}\mathrm{f}}=\sqrt{{Z}_{\mathrm{2}}^{\mathrm{2}}-{Z}_{\mathrm{2}\mathrm{p}}^{\mathrm{2}}}$ where *Z*_{2p} is the contribution of the
particle-scale roughness. Using the data described in
Sect. 4.1, values of *Z*_{2p}≈0.23–0.24 are
obtained. This would result in a correction of the mean JRC values for the
different groups of surfaces shown in Fig. 6a from ≈23.7 to ≈22.1 for the smallest and from ≈32.2 to ≈31.8 for the largest values of the JRC. This shows that the potential
corrections are not significant and, in most cases, well inside the scatter of
the calculated JRC values. In addition, we did run two small sets of simulations
using a wider range of particle sizes than the “standard” models described
in Sect. 3.1, i.e., a larger ratio between maximum and minimum
particle radius (${R}_{\text{max}}:{R}_{\text{min}}=\mathrm{1.0}:\mathrm{0.15}$ and
${R}_{\text{max}}:{R}_{\text{min}}=\mathrm{1.0}:\mathrm{0.1}$), to see if the particle size range
had any effect on the surface properties. These sets consisted of five
simulations each, all performed under true triaxial conditions using
*σ*_{2}=6 MPa and *σ*_{3}=0. The results did not show a
statistically significant difference in Hurst exponent or JRC compared to the
equivalent simulations performed using the particle radius range
${R}_{\text{max}}:{R}_{\text{min}}=\mathrm{1.0}:\mathrm{0.2}$ (Table 1).

In numerical models, there is a slightly higher anisotropy in the models with
transversely isotropic confinement (*σ*_{2}=*σ*_{3}) of up to
8 % difference in JRC between the directions, whereas in the models
with *σ*_{2}≠*σ*_{3} the difference is less than 3 % in all
cases. In the rock samples, which are also deformed under conditions where
*σ*_{2}=*σ*_{3}, the anisotropy is much higher; i.e., the ratio between
the JRC in the two directions is $\approx \mathrm{1}:\mathrm{1.2}$ in the sandstone and $\approx \mathrm{1}:\mathrm{2.3}$ in the limestone. However, due to small number of fracture surfaces
available for analysis from the laboratory experiments, it is not clear if this
strong anisotropy, and the large difference between the limestone and the
sandstone sample, is a general property of fracture surfaces generated under
comparable conditions or just an artifact of the specific samples studied. In
general, the strong anisotropy which was observed in the laboratory
experiments, in particular in the limestone, was not replicated in the
numerical models. The reason for the stronger directional anisotropy in the
natural rock samples is not clear yet. A key difference between the
microscale mechanics of laboratory and numerical experiments is that the
natural rocks can undergo grain size reduction during the fracture process,
whereas this mechanism is not implemented in the numerical models used in this
paper. This might explain why the numerical models, at least in our
experiments, do not produce the striations observed in the natural rock
samples. A possibility to test this hypothesis in future work would be to
extend the numerical models to use breakable particle clusters to represent
rock grains instead of single particles. This approach has been shown to yield
insights into the micromechanics of grain size reduction processes, for
example, in fault gouge (Abe and Mair, 2005; Mair and Abe, 2008, 2011) and in
compression experiments (Thornton et al., 2004). However, it also
significantly increases the required computational effort for the
simulations. A computationally less expensive option to include grain size
reduction into the numerical models might be to adapt the empirical particle
replacement approach developed by Cleary (2001) to the specific
requirements of the simulation of rock fracture under triaxial
loading. However, as Weerasekara et al. (2013) point out, this approach is
strongly dependent on the availability of good calibration data for the grain
fracture under the specific stress and strain rate conditions of the process
modeled. Further insights could also be provided by additional laboratory
experiments, for example, to test if the difference in anisotropy between
numerical and laboratory experiments also exists under true triaxial
conditions (*σ*_{2}≠*σ*_{3}).

Smoothing due to abrasion while sliding is, in general, an important mechanism
for the modification of rough surfaces. In particular, slip along the surface
can result in a significant reduction of the Hurst exponent for profiles
parallel to the slip direction down to values below 0.5
(Candela et al., 2012, Table 1b). However, those large reductions appear to
apply mainly to faults with large amounts of slip, i.e., several meters up to
kilometers. In contrast, data from laboratory experiments published in the
literature (Amitrano and Schmittbuhl, 2002; Davidesko et al., 2014; Badt et al., 2016) suggest that this process is unlikely to have a sufficiently
large effect at the small shear offsets in both numerical models and
experimental samples studied here to explain the observed differences. To
investigate if the roughness evolution of the fracture surfaces with
increasing deformation of the sample plays a role in our numerical model, we
did perform a small number of simulations which did not stop immediately after
the formation of the fractures but instead continued deformation to a total
axial strain of up to 12 %. This is significantly larger than the
strain occurring in the laboratory experiments, where total axial shortening
did not exceed about 2 %. In particular, the amount of shortening
occurring after the peak axial stress was reached, i.e., after failure, was
generally less than 1 %. The obtained Hurst exponents did show no
significant trend with increasing strain of the model and offset of the shear
fracture (Fig. S4 in the Supplement). While the average of the Hurst exponents
from the six surfaces investigated could be considered as showing a slight
increasing trend for axial strains up to 8 % (Fig. S5 in the
Supplement), the increase of 0.03 is about an order of magnitude too small to
explain the observed differences between numerical and experimental
surfaces. However, it would be compatible with the effect observed by
Amitrano and Schmittbuhl (2002). For one of the models, we also calculated
the JRC of the surfaces at various stages of the simulation. The data show
that there is also no significant change of the JRC for the shear offset
considered in this model, which would be equivalent to ≈1 cm
in the laboratory samples, and under the conditions of this model, i.e., true
triaxial stress with *σ*_{2}=7.5 MPa, *σ*_{3}=3 MPa
(Fig. S6 in the Supplement). This seems to confirm again that under the small
shear offsets relevant for our experiments, there is very little evolution of
the surface roughness, at least as far as it concerns the roughness parameters
calculated here (Hurst exponent, JRC). In particular, the data would suggest
that any effects due to the small, but the non-zero, shear offset in the
laboratory experiments is much too small to explain the observed differences
between numerical simulations and laboratory experiments.

Based on the results from the numerical models, there appears to be a trend
towards higher roughness for fracture surfaces generated under transversely
isotropic stress conditions, i.e., standard triaxial compression (${\mathit{\sigma}}_{\mathrm{1}}>{\mathit{\sigma}}_{\mathrm{2}}={\mathit{\sigma}}_{\mathrm{3}}$) compared to those generated under true triaxial
conditions (*σ*_{2}≠*σ*_{3}). This trend was shown for both
geometrical roughness measures used in the analysis of the data from the
numerical experiments, i.e., the JRC
(Fig. 6) and also the rms roughness
(Fig. 8). A possible, but at this stage purely speculative,
idea to explain this observation might be that, if we assume that the
through-going fractures, which we analyze, form by coalescence from smaller,
precursory fractures, those precursory fractures have their strike angles
constrained to a narrow range if *σ*_{2}≠*σ*_{3}, but that there is
no such constraint if *σ*_{2}=*σ*_{3}. If this is the case, then the
coalescence of those precursory fractures might lead to smoother large-scale
surfaces if they all have similar orientations compared to when they have
random strike directions. Unfortunately, the numerical models used in this work
do not have the resolution necessary to test this hypothesis.

Additionally, a difference in the roughness between the surfaces on tensile and compressive (i.e., shear) fractures generated under unconfined conditions has been observed, with the tensile fractures showing a smaller roughness. This effect appears to be more pronounced if the roughness is measured in terms of the JRC compared to the rms roughness. Should these effects be confirmed by further work, and in particular by comparison with more experimental data, it could be used to provide additional input data to, for example, permeability estimations of fracture networks or geomechanical fault stability calculations.

The analysis of the roughness scaling properties of the surfaces in terms of
the height–height correlation function shows that the fracture surfaces
generated in the numerical models are self-affine with Hurst exponents around
0.3–0.45. This value is in disagreement with the majority of field and
experimental studies (Bouchaud et al., 1990; Schmittbuhl et al., 1993, 1995; Bouchaud, 1997) which find a “universal” Hurst
exponent *H*≈0.8. However, low Hurst exponents in the range *H*≈0.4–0.5 have previously also been found in other numerical models of the
generation of rough fractures such as 3-D random fuse networks
(Alava et al., 2006).

The Hurst exponents of the surfaces generated in the numerical models can be corrected for the influence of the particle-scale roughness in a similar way to the procedure described above for the correction of the joint roughness coefficients. It would require correcting the rms roughness values in the height–height correlation function for each individual distance bin and obtaining a power-law fit based on the corrected data points (Fig. 16). However, while these corrections do lead to slightly higher calculated Hurst exponents, the increase is at most about 0.05 and therefore the effect is far too small to explain the discrepancy.

The data obtained from the fracture surfaces generated in triaxial tests on
the limestone sample (*H*≈0.75) are compatible with this “universal
exponent”. In contrast, the sandstone sample shows a lower Hurst exponent
(*H*≈0.6) than the limestone sample but not as low as the numerical
models. There are experimental data for sandstone in the literature showing
Hurst exponents even lower than our sandstone sample and in fact close to the
results from the numerical models, i.e., $H=\mathrm{0.47}\pm \mathrm{0.04}$,
(Boffa et al., 1998) and similar data from a synthetic, sandstone-like
material made from sintered glass beads ($H=\mathrm{0.40}\pm \mathrm{0.04}$,
Ponson et al., 2006). Both those studies investigated tensile (mode-1)
fractures. Boffa et al. (1998) used a direct tension setup with a
pre-notched sample to initiate the fracture at a defined location, whereas
Ponson et al. (2006) used a modified Brazilian test where a compressive
load is applied to two opposite points on the circumference of the cylindrical
sample to generate a tensile stress in the stress in the central part of the
disk (Jaeger et al., 2007; Fjaer et al., 2008). However, our numerical models
do not show a dependence of the Hurst exponent on the fracture mode
(Fig. 11).

Nigon et al. did observe a transition from a Hurst exponent of 0.74 to a lower value of 0.5 below a length scale of about 0.1 mm in natural joint surfaces in sandstone (Nigon et al., 2017, Fig. 9). However, this transition scale from a “jointing induced roughness” to a “grain induced roughness” is at a scale comparable to the mean grain size in their material. The equivalent length scale in our numerical models would be the mean particle diameter, i.e., below 1 model unit, which is well below the length range used to fit the scaling law (Fig. 9). This difference in scales shows that the Hurst-exponents in our numerical models are completely calculated above the “transition scale” of Nigon et al. (2017) and therefore should belong to the regime described as “jointing induced roughness” by them. This means that the low values of the Hurst-exponents in the numerical cannot be explained by the “grain induced roughness” regime of Nigon et al. (2017).

When comparing the data from the numerical models to the relation between
fractal dimension *D* and JRC proposed by Ficker (2017),
i.e., $\text{JRC}\approx \mathrm{50}(D-\mathrm{1})$, the surfaces show on average a slightly
smaller JRC than would be expected based on their fractal dimension *D*
calculated from the Hurst exponent as $D=\mathrm{2}-H$
(Fig. 17). Interestingly, the data from the sandstone sample
plot even further below the relation by Ficker (2017). The data from
the limestone sample are difficult to compare due to the large anisotropy of
the JRC and are therefore not plotted in Fig. 17.

It has been suggested by Ponson et al. (2007) that the observed Hurst
exponent is an indicator for the failure mode, *H*≈0.8 for “damage
fracture”, i.e., coalescence from microcracks and *H*≈0.4 for
“brittle fracture”, i.e., continuous propagation of the crack. However, we
have not been able to confirm this for our numerical experiments. Looking at
the relative timing of bonds breaking suggests that the fracture surfaces in
the DEM models grow by coalescence of microcracks despite having a Hurst
exponent closer to 0.4. For examples of the general evolution of the
microcrack distribution, see Figs. S7 and S8 in the Supplement.

The dependence of the variability of the measured Hurst exponent on the size of the analyzed surface on both limestone and sandstone samples suggests the large scatter observed in the Hurst exponents from the numerical models could be a resolution issue. The sandstone sample has a maximum grain size of about 200 µm. This results in a ratio between the length and width of the analyzed fracture surface and the maximum grain size of between 250:1 and 350:1, whereas this ratio is only in the range between 30:1 and 60:1 in the numerical models. The limestone sample is even more fine grained than the sandstone sample.

Amitrano and Schmittbuhl (2002) find a weak decrease of the roughness
exponent with increasing confinement if no further shear displacement is
imposed on the surfaces after fracture. This is similar to the trend observed
in our numerical simulation data (Fig. 11), although at
different absolute values of the Hurst exponent, which are in the range
between 0.3 to 0.45 in our data and between 0.7 to 0.77 in
Amitrano and Schmittbuhl (2002). Also, this stress dependence cannot be
directly compared because of differences in the mechanical properties between
the simulated material in our case and the real granite.
Amitrano and Schmittbuhl (2002) do not explicitly give the unconfined
compressive strength (UCS) of the granite. Extrapolating from their Fig. 3
suggests a value of around 300 MPa, although a calculation from their
internal cohesion (37 MPa) and friction angle (55±2^{∘})
gives a value closer to 240 MPa. Combined with the confining stress
used in their work of *σ*_{3}≈20–80 MPa, this suggests that
the ratio between UCS and the confining stress is in a similar range as in the
numerical models used here, where UCS=80 MPa and *σ*_{3}=0–15 MPa.

Synthetic fracture surfaces have been generated in numerical simulations of rock deformation experiments using the discrete element method (DEM). Results of a statistical analysis demonstrate that the generated surfaces are self-affine. Further analysis has shown no dependency of roughness measures such as rms roughness and the joint roughness coefficient (JRC) on the confining stress. One exception is the observation that samples fractured under true anisotropic conditions (${\mathit{\sigma}}_{\mathrm{1}}>{\mathit{\sigma}}_{\mathrm{2}}>{\mathit{\sigma}}_{\mathrm{3}}$) show lower JRC and lower rms roughness than samples fractured under transversal isotropic confinement (${\mathit{\sigma}}_{\mathrm{1}}>{\mathit{\sigma}}_{\mathrm{2}}={\mathit{\sigma}}_{\mathrm{3}}$), at least for numerical models. For natural rock samples this effect has not been tested yet. Photogrammetric analysis of shear fracture surfaces on two rock samples has shown that the choice of sampling area can influence the roughness data obtained. Results show, for example, a variation of ±0.1 in the Hurst exponent between small sampling areas on the same surface of a rock sample.

Comparing the numerical results with laboratory experiments and additional data obtained from the literature suggests that the trends observed in the numerical parameter study are valid, but it also shows some discrepancies in the absolute values of some of the roughness parameters. In particular, the fracture surfaces generated in the DEM simulations show a higher joint roughness coefficient compared to natural rock samples and a lower Hurst exponent. The comparison also shows a stronger directional anisotropy of the roughness in the real rock samples compared to the numerical simulations. The reason for this result is not clear so far and should be subject to further investigation. One possible cause might be the occurrence of grain size reduction in real rocks, which is not implemented in the current numerical models.

The open-source code for DEM simulations is available at https://launchpad.net/esys-particle (last access: 13 October 2021).

The scripts to reproduce the input data sets are available upon request from the authors.

The supplement related to this article is available online at: https://doi.org/10.5194/se-12-2407-2021-supplement.

SA performed the numerical simulations and the data analysis on the numerical models, and wrote the initial draft of the manuscript; HD performed the analysis of the laboratory samples and edited the manuscript.

The contact author has declared that neither they nor their co-author have any competing interests.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

The triaxial deformation tests were carried out at Technische Universität Darmstadt. The authors thank Jessica McBeck, Francesco Salvini and an anonymous referee for reviewing this paper.

The work was carried out within the project “PERMEA”, funded by the German Federal Ministry of Education and Research (BMBF) (funding ID FKZ 03G0865B).

This paper was edited by Federico Rossetti and reviewed by Francesco Salvini, Jessica McBeck, and one anonymous referee.

Abe, S. and Mair, K.: Grain fracture in 3D numerical simulations of granular shear, Geophys. Res. Lett., 32, L05305, https://doi.org/10.1029/2004GL022123, 2005. a, b

Abe, S., Place, D., and Mora, P.: A Parallel Implementation of the Lattice Solid Model for the Simulation of Rock Mechanics and Earthquake Dynamics, Pure Appl. Geophys., 161, 2265–2277, https://doi.org/10.1007/s00024-004-2562-x, 2004. a

Abe, S., van Gent, H., and Urai, J. L.: DEM simulation of normal faults in cohesive materials, Tectonopysics, 512, 12–21, https://doi.org/10.1016/j.tecto.2011.09.008, 2011. a

Ahmadi, M., Taleghani, A. D., and Sayers, C.: The effects of roughness and offset on fracture compliance ratio, Geophys. J. Int., 205, 454–463, https://doi.org/10.1093/gji/ggw034, 2016. a

Alava, M. J., Nukala, P. K. V. V., and Zapperi, S.: Statistical models of fracture, Adv. Phys., 55, 349–476, https://doi.org/10.1080/00018730300741518, 2006. a, b

Amanatides, J. and Woo, A.: A Fast Voxel Traversal Algorithm for Ray Tracing, in: Proceedings of EuroGraphics, Eurographics Association, Geneva, 87, 1987. a

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

Angheluta, L., Candela, T., Mathiesen, J., and Renard, F.: Effect of Surface Morphology on the Dissipation During Shear and Slip Along a Rock–Rock Interface that Contains a Visco-elastic Core, Pure Appl. Geophys., 168, 2335–2344, https://doi.org/10.1007/s00024-011-0272-8, 2011. a

Badt, N., Hatzor, Y. H., Toussaint, R., and Sagy, A.: Geometrical evolution of interlocked rough slip surfaces: The role of normal stress, Earth Planet. Sc. Lett., 443, 153–156, https://doi.org/10.1016/j.epsl.2016.03.026, 2016. a

Barton, N.: Review of a new shear-strength Criterion for Rock Joints, Eng. Geol., 7, 287–332, 1973. a, b, c, d

Barton, N. and Choubey, V.: The Shear Strength of Rock Joints in Theory and Practice, Rock Mech., 10, 1–54, 1977. a, b, c, d, e

Bisdom, K., Bertotti, G., and Nick, H. M.: A geometrically based method for predicting stress-induced fracture aperture and flow in discrete fracture networks, AAPG Bull., 100, 1075–1097, https://doi.org/10.1306/02111615127, 2016. a

Boffa, J., Allain, C., and Hulin, J.: Experimental analysis of fracture rugosity in granular and compact rocks, Eur. Phys. J.-Appl. Phys., 2, 281–289, https://doi.org/10.1051/epjap:1998194, 1998. a, b

Bouchaud, E.: Scaling properties of cracks, J. Phys.-Condens. Mat., 9, 4319–4344, https://doi.org/10.1088/0953-8984/9/21/002, 1997. a, b, c, d

Bouchaud, E., Lapasset, G., and Planes, J.: Fractal Dimension of Fractured Surfaces: a Universal Value, Europhys. Lett., 13, 73–79, 1990. a

Briggs, S., Karney, B. W., and Sleep, B. E.: Numerical modeling of the effects of roughness on flow and eddy formation in fractures, Int J. Rock Mech. Geotech. Eng., 9, 105–115, https://doi.org/10.1016/j.jrmge.2016.08.004, 2017. a

Brown, S. R. and Scholz, C. H.: Broad bandwidth study of the topography of natural rock surfaces, J. Geophys. Res., 90, 12575–12582, https://doi.org/10.1029/jb090ib14p12575, 1985. 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, https://doi.org/10.1007/s00024-009-0521-2, 2009. a, b, c, d, e, f

Candela, T., Renard, F., Bouchon, M., Schmittbuhl, J., and Brodsky, E. E.: Stress Drop during Earthquakes: Effect of Fault Roughness Scaling, B. Seismol. Soc. Am., 101, 2369–2387, https://doi.org/10.1785/0120100298, 2011a. a

Candela, T., Renard, F., Schmittbuhl, J., Bouchon, M., and Brodsky, E. E.: Fault slip distribution and fault roughness, Geophys. J. Int., 197, 959–968, https://doi.org/10.1111/j.1365-246X.2011.05189.x, 2011b. 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., 117, B08409, https://doi.org/10.1029/2011JB009041, 2012. a, b, c

Chen, Z., Narayan, S., Yang, Z., and Rahman, S.: An experimental investigation of hydraulic behaviour of fractures and joints in granitic rock, Int. J. Rock Mech. Min., 37, 1061–1071, 2000. a

Cleary, P.: Recent Advances in DEM Modelling of Tumbling Mills, Miner. Eng., 14, 1295–1319, 2001. a

Cundall, P. A. and Strack, O.: A discrete numerical model for granular assemblies, Géotechnique, 29, 47–65, 1979. a, b, c

Davidesko, G., Sagy, A., and Hatzor, Y. H.: Evolution of slip surface roughness through shear, Geophys. Res. Lett., 41, 1492–1498, https://doi.org/10.1002/2013GL058913, 2014. a

Donze, F., Mora, P., and Magnier, S.-A.: Numerical simulation of faults and shear zones, Geophys. J. Int., 116, 46–52, 1994. a, b

El-Soudani, S. M.: Profilometric analysis of fracture surfaces, Metallography, 11, 247–336, https://doi.org/10.1016/0026-0800(78)90045-9, 1978. a

Fakhimi, A. and Gharahbagh, E. A.: Discrete element analysis of the effect of pore size and pore distribution on the mechanical behavior of rock, Int. J. Rock Mech. Min.., 48, 77–85, https://doi.org/10.1016/j.ijrmms.2010.08.007, 2011. a

Faulkner, D., Jackson, C., Lunn, R., Schlische, R., Shipton, Z., Wibberley, C., and Withjack, M.: A review of recent developments concerning the structure, mechanics and fluid flow properties of fault zones, J. Struct. Geol., 21, 1557–1575, https://doi.org/10.1016/j.jsg.2010.06.009, 2010. a

Ficker, T.: Fractal properties of joint roughness coefficients, Int. J. Rock Mech. Min., 94, 27–31, https://doi.org/10.1016/j.ijrmms.2017.02.014, 2017. a, b, c

Fjaer, E., Holt, R., Horsrud, P., Raaen, A., and Risnes, R.: Petroleum related rock mechanics, 2nd Edn., Elsevier, Amsterdam, 2008. a

Griffith, W. A., Nielsen, S., Toro, G. D., and Smith, S. A. F.: Rough faults, distributed weakening, and off-fault deformation, J. Geophys. Res., 115, B08409, https://doi.org/10.1029/2009JB006925, 2010. a

Hobbs, B.: The significance of structural geology in rock mechanics, in: Comprehensive Rock Engineering, Pergamon Press, Oxford, Vol. 1, 25–62, 1993. a

Jaeger, J., Cook, N., and Zimmerman, R.: Fundamentals of rock mechanics, 4th Edn., Blackwell, Oxford, UK, 2007. 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, https://doi.org/10.1016/j.ijheatmasstransfer.2016.10.010, 2017. a

Kottwitz, M. O., Popov, A. A., Baumann, T. S., and Kaus, B. J. P.: The hydraulic efficiency of single fractures: correcting the cubic law parameterization for self-affine surface roughness and fracture closure, Solid Earth, 11, 947–957, https://doi.org/10.5194/se-11-947-2020, 2020. a, b

Kottwitz, O.: Scale Invariant Roughness Quantification and Anisotropy Exposure of Rock Discontinuities, Master's thesis, Johannes-Gutenberg-Universität Mainz, 2017. a

Koyama, T. and Jing, L.: Effects of model scale and particle size on micro-mechanical properties and failure processes of rocks – A particle mechanics approach, Eng. Anal. Bound. Elem., 31, 458–472, https://doi.org/10.1016/j.enganabound.2006.11.009, 2007. a

Li, Y. and Zhang, Y.: Quantitative estimation of joint roughness coefficient using statistical parameters, Int. J. Rock Mech. Min., 77, 27–35, https://doi.org/10.1016/j.ijrmms.2015.03.016, 2015. a, b, c, d, e, f, g, h

Mair, K. and Abe, S.: 3D numerical simulations of fault gouge evolution during shear: Grain size reduction and strain localization, Earth Planet. Sc. Lett., 274, 72–81, 2008. a, b

Mair, K. and Abe, S.: Breaking Up: Comminution Mechanisms in Sheared Simulated Fault Gouge, Pure Appl. Geophys., 168, 2277–2288, https://doi.org/10.1007/s00024-011-0266-6, 2011. a

Mandelbrot, B. B.: Self-Affine Fractals and Fractal Dimension, Phys. Scripta, 32, 257–260, https://doi.org/10.1088/0031-8949/32/4/001, 1985. a, b, c

Mandelbrot, B. B. and Van Ness, J. W.: Fractional Brownian Motions, Fractional Noises and Applications, SIAM Rev., 10, 422–437, https://doi.org/10.1137/1010093, 1968. a

Mora, P. and Place, D.: Simulation of the Frictional Stick-slip Instability, Pure Appl. Geophys., 143, 61–87, 1994. a, b

Myers, N. O.: Charackterization of Surface Roughness, Wear, 5, 182–189, 1962. a

Nigon, B., Englert, A., Pascal, C., and Saintot, A.: Multiscale Characterisation of Joint Surface Roughness, J. Geophys. Res., 122, 9714–9728, https://doi.org/10.1002/2017JB014322, 2017. a, b, c

Ohagen, M.: Untersuchungen zur Rauheit von Klüften in DEM-Modellen, BSc-Thesis, Johannes-Gutenberg-Universität Mainz, 2019. a, b, c, d

Okubo, P. G. and Dietrich, J. H.: Effects of physical fault properties on frictional instabilities produced on simulated faults, J. Geophys. Res., 89, 5817–5827, 1984. a

Place, D. and Mora, P.: A random lattice solid model for simulation of fault zone dynamics and fracture process, in: Bifurcation and Localisation Theory for Soils and Rocks 99, edited by: Mühlhaus H.-B., Dyskin, A. V., and Pasternak, E., AA Balkema Rotterdam/Brookfield, 2001. a

Ponson, L., Auradou, H., Vié, P., and Hulin, J.-P.: Low Self-Affine Exponents of Fractured Glass Ceramics Surfaces, Phys. Rev. Lett., 97, 125501, https://doi.org/10.1103/PhysRevLett.97.125501, 2006. a, b

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, 03618, https://doi.org/10.1103/PhysRevE.76.036108, 2007. a, b

Power, W. L. and Durham, W. B.: Topography of Natural and Artificial Fractures in Granitic Rocks: Implications for Studies of Rock Friction and Fluid Migration, Int. J. Rock Mech. Min., 34, 979–989, 1997. a

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, https://doi.org/10.1029/2005GL025038, 2006. a

Sayles, R. S. and Thomas, T. R.: The Spatial Representation of Surface Roughness by Means of the Structure Function: A practical Alternative to Correlation, Wear, 42, 263–276, 1977. a

Schmittbuhl, J., Gentier, S., and Roux, S.: Field measurements of the roughness of fault surfaces, Geophys. Res. Lett., 20, 639–641, https://doi.org/10.1029/93gl00170, 1993. a, b

Schmittbuhl, J., Schmitt, F., and Scholz, C.: Scaling invariance of crack surfaces, J. Geophys. Res., 100, 5953–5973, https://doi.org/10.1029/94jb02885, 1995. a, b

Schöpfer, M. P. J., Abe, S., Childs, C., and Walsh, J. J.: The impact of porosity and crack density on the elasticity, strength and friction of cohesive granular materials: Insights from DEM modelling, Int. J. Rock Mech. Min., 46, 250–261, 2009. a

Schöpfer, M. P. J., Arslan, A., Walsh, J. J., and Childs, C.: Reconciliation of contrasting theories for fracture spacing in layered rocks, J. Struct. Geol., 33, 551–565, https://doi.org/10.1016/j.jsg.2011.01.008, 2011. a

Thornton, C., Ciomocos, M., and Adams, M.: Numerical simulations diametrical compression tests on agglomerates, Powder Technol., 140, 258–267, https://doi.org/10.1016/j.powtec.2004.01.022, 2004. a

Tse, R. and Cruden, D. M.: Estimating Joint Roughness Coefficients, Int. J. Rock. Mech. Min., 16, 303–307, 1979. a, b, c

Turcotte, D. L.: Fractals, chaos, self-organized citicality and tectonics, TerraNova, 4, 4–12, 1992. a

Wang, Y., Abe, S., Latham, S., and Mora, P.: Implementation of particle-scale rotation in the 3-D lattice solid model, Pure Appl. Geophys., 163, 1769–1785, 2006. a

Watanabe, N., Hirano, N., and Tsuchiya, N.: Determination of aperture structure and fluid flow in a rock fracture by high-resolution numerical modeling on the basis of a flow-through experiment under confining pressure, Water Resour. Res., 44, W06412, https://doi.org/10.1029/2006WR005411, 2008. a

Weatherley, D.: Investigations on the role of microstructure in brittle failure using discrete element simulations, Geophysical Research Abstracts, 13, EGU2011-9476, 2011. a

Weerasekara, N., Powell, M., Cleary, P., Tavares, L., Evertsson, M., Morrison, R., Quist, J., and Carvalho, R.: The contribution of DEM to the science of comminution, Powder Technol., 248, 3–24, https://doi.org/10.1016/j.powtec.2013.05.032, 2013. a

Yang, Z. and Lo, S.: An Index for Describing the Anisotropy of Joint Surfaces, Int. J. Rock Mech. Min., 34, 1031–1044, https://doi.org/10.1016/S1365-1609(97)80012-3, 1997. a

Yoon, J. S., Zang, A., and Stephansson, O.: Simulating fracture and friction of Aue granite under confined asymmetric compressive test using clumped particle model, Int. J. Rock Mech. Min., 49, 68–83, https://doi.org/10.1016/j.ijrmms.2011.11.004, 2012. a

Yu, X. and Vayssade, B.: Joint Profiles and their Roughness Parameters, Int. J. Rock Mech. Min., 28, 333–336, 1991. a, b

Zambrano, M., Pitts, A. D., Salama, A., Volatili, T., Giorgioni, M., and Tondi, E.: Analysis of Fracture Roughness Control on Permeability Using SfM and Fluid Flow Simulations: Implications for Carbonate Reservoir Characterization, Geofluids, 2019, ID 4132386, https://doi.org/10.1155/2019/4132386, 2019. a

Zoback, M. D.: Reservoir Geomechanics, Cambridge University Press, Cambridge, UK, 2007. a