Roughness of fracture surfaces in numerical models and laboratory experiments

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.

1 Introduction

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₂ (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₂, 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 Method

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).

Figure 1 Numerical modeling workflow. (a) DEM specimen used for deformation experiments. Colors show particle size, purple arrows symbolize confining stress, gray arrows show compression direction. (b) Fragment identification in fractured DEM model. Colors show fragment size (red – large, blue – small). Red parts (top left and bottom right) show two major fragments, blue/white (i.e., fine-grained) material along diagonal shows shear zone. (c) Two major fragments extracted from fractured DEM model. Colors show fragment size (volume). (d) Fragment extracted from DEM model. Rough fracture surface visible. (e) Point cloud surface generated from DEM model. Outer surfaces of the initial DEM specimen visible right and bottom. (f) Filtered point cloud used for analysis. Non-fracture surfaces and outlying points removed.

Figure 2 Simplified 2-D sketch of the ray-casting method. The gray particles are assumed to belong to the same fragment of the deformed sample, and the black and gray crosses show the fragment surface calculated from the line–particle intersections using multiple view directions. Black lines and the black arrow show primary view direction; light gray lines and arrows show additional view directions at a 30^∘ angle to the primary direction.

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

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }(x_{\mathrm{0}},y_{\mathrm{0}},z_{\mathrm{0}})={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{0}}^n(x_i,y_i,z_i),\end{array}$$

and its orientation is determined by the two major eigenvectors e₁ and e₂ of the covariance matrix C of the point cloud. The third eigenvector of the covariance matrix then determines the normal n_fp=e₃ 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

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }{x}^{\prime }=(\mathit{p}-\mathit{b})\cdot {\mathit{e}}_{\mathrm{1}}\text{(3)}& {\displaystyle }& {\displaystyle }{y}^{\prime }=(\mathit{p}-\mathit{b})\cdot {\mathit{e}}_{\mathrm{2}}\text{(4)}& {\displaystyle }& {\displaystyle }{z}^{\prime }=(\mathit{p}-\mathit{b})\cdot {\mathit{e}}_{\mathrm{3}}.\end{array}$$

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).

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }{\text{JRC}}_{\mathrm{1}}=\mathrm{32.2}+\mathrm{32.47}\log \left(Z_{\mathrm{2}}\right)\text{(6)}& {\displaystyle }& {\displaystyle }{\text{JRC}}_{\mathrm{31}}=\mathrm{558.68}\sqrt{R_{\mathrm{p}}}-\mathrm{557.13}\text{(7)}& {\displaystyle }& {\displaystyle }{\text{JRC}}_{\mathrm{34}}=\mathrm{92.97}\sqrt{\mathit{\delta }}-\mathrm{5.25},\end{array}$$

where R_p is the “roughness profile index”, $\mathit{\delta }=R_{\mathrm{p}}-\mathrm{1}$ the “profile elongation index” and Z₂ 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

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }{R}_{\mathrm{p}}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^{N-\mathrm{1}}\sqrt{(x_{i+\mathrm{1}}-x_i{)}^{\mathrm{2}}+(y_{i+\mathrm{1}}-y_i{)}^{\mathrm{2}}}}{{\sum }_{i=\mathrm{1}}^{N-\mathrm{1}}\sqrt{x_{i+\mathrm{1}}-x_i}}}\end{array}$$

and

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }{Z}_{\mathrm{2}}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^{N-\mathrm{1}}{\displaystyle \frac{(y_{i+\mathrm{1}}-y_i{)}^{\mathrm{2}}}{x_{i+\mathrm{1}}-x_i}}},\end{array}$$

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

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }(x,y,z)\to (ax,ay,a^{H}z),\end{array}$$

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).

Figure 3 Height and distance relations of points in the point cloud used to calculate the height–height correlation function. (a) Arrangement of points above a fitted mean plane: dashed grid showing mean plane, black lines symbolizing orthogonal distance between plane and points, and red/green lines showing relative orientation between points. (b) Distance (Δr) and height difference (Δz) between points.

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

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }c\left(\mathrm{\Delta }r\right)=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum _{\mathrm{\Delta }r-w/\mathrm{2}}^{\mathrm{\Delta }r+w/\mathrm{2}}\mathrm{\Delta }{z}^{\mathrm{2}}},\end{array}$$

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 }=\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.

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }c(\mathrm{\Delta }r,\mathit{\varphi })=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum _{\begin{array}{c}\mathrm{\Delta }r-w_r/\mathrm{2}\\ \mathit{\varphi }-w_{\mathit{\varphi }}/\mathrm{2}\end{array}}^{\begin{array}{c}\mathrm{\Delta }r+w_r/\mathrm{2}\\ \mathit{\varphi }+w_{\mathit{\varphi }}/\mathrm{2}\end{array}}\mathrm{\Delta }{z}^{\mathrm{2}}},\end{array}$$

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.

3 Experiments

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.

Figure 4 Types of fracture surfaces studied in this work: (a) numerical (DEM) model, (b) limestone fragment generated in triaxial deformation experiment, (c) sandstone fragment generated in uniaxial deformation experiment.

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 (σ₁>0, ${\mathit{\sigma }}_{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{3}}=\mathrm{0}$), unconfined tension (σ₃<0, ${\mathit{\sigma }}_{\mathrm{1}}={\mathit{\sigma }}_{\mathrm{2}}=\mathrm{0}$), standard triaxial compression (σ₁>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, σ₁ is parallel to the y axis, σ₂ is parallel to the x axis, and σ₃ is parallel to the z axis, whereas in the unconfined tensile models σ₃, i.e., the extension direction, is parallel to the y axis, σ₁ is parallel to the x axis, and σ₂ 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=σ₂ and σ_zz=σ₃. 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 σ₁=σ₂. 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⁻³.

In all deformation experiments, the final loading plate speed was set to $\approx \mathrm{17}\mathrm{cm}{\mathrm{s}}^{-\mathrm{1}}$. This is significantly higher than in real experiments, but using real lab values ($\mathrm{\mu }\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}\mathrm{\dots }\mathrm{mm}{\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.

Figure 5 Confining stress range covered by the numerical models, combining the experiments in this work and the data from Ohagen (2019). (a) Numerical models only, using absolute stress values. (b) Numerical and laboratory experiments, using stress values scaled by the unconfined compressive strength of the respective material. Hatched segment in top left of the diagrams: parameter space excluded by the condition σ₃≤σ₂.

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⁻² 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 Results

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 σ₁, i.e., 55–65^∘ assuming σ₁ to be vertical. The strike direction of the majority of the fractures was within ≈10^∘ of σ₂ in the true triaxial models (σ₂≠σ₃) and more or less randomly distributed under transverse isotropic stress conditions (σ₂=σ₃).

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₂ (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₂≈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 (σ₂=σ₃) generally have higher JRC values than models fractured under true triaxial conditions, i.e., σ₂≠σ₃. The directly calculated geometric roughness measures, i.e., R_p and Z₂, show a very similar pattern (Fig. 6b and c).

Figure 6 Geometric roughness measures for surfaces generated at different stress conditions in DEM models. Black: True triaxial compression, bed: Standard triaxial compression (transverse isotropic confinement), blue: unconfined extension. (a) Approximated JRC values calculated based on Eq. (5). (b) “Root mean square of the first deviation” Z₂ (Eq. 9). (c) “Profile elongation index” R_p (Eq. 8). Error bars show standard deviation.

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.

Figure 7 Distribution of heights of a simulated fracture surface above a “best-fit” plane. Data are taken from a model with σ₂=5 MPa and σ₃=0 MPa. (a) Map view of the surface colored by height above the “best-fit” plane. (b) Probability density of heights and fitted normal distribution.

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 (σ₂=σ₃) and true triaxial (σ₂≠σ₃) 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.

Figure 8 Average rms roughness values for surfaces generated at different stress conditions: black: true triaxial compression, red: standard triaxial compression (transverse isotropic confinement), blue: unconfined extension. Error bars show standard deviation.

Figure 9 Log–log plot of the height–height correlation function of the two surfaces of a single fracture. Red and black symbols show the rms height differences calculated for each distance bin for the two surfaces. The straight lines are fitted to the linear section of the data in log–log space, showing a power-law dependence.

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.

Figure 10 Log–log plot of the height–height correlation function of a fracture surface generated in a numerical deformation experiment (red diamonds) and a “quasi-planar” surface generated by cutting the particle packing used in the experiment with a plane (black circles).

Figure 11 Average and variability of Hurst exponents for surfaces generated with different average confining stress $({\mathit{\sigma }}_{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{3}})/\mathrm{2}$. Black: triaxial compression models, this work, blue: unconfined extension models, this work, red and pink: triaxial compression models; data are from Ohagen (2019). Error bars show top and bottom quartiles.

Figure 12 Average and variability of Hurst exponents for surfaces generated with different ratio of confining stresses σ₂ and σ₃. Black: this work, red: data from Ohagen (2019). Error bars show top and bottom quartiles.

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² for the sandstone and ≈25 cm² 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.

Figure 13 Distribution of heights for fracture surfaces generated in laboratory deformation experiments. Top: limestone sample deformed at ${\mathit{\sigma }}_{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{3}}=\mathrm{30}\mathrm{MPa}$, bottom: sandstone sample deformed in unconfined compression test. (a) Map view of the limestone surface colored by height above the “best-fit” plane. (b) Probability density of heights and fitted normal distribution for limestone surface (c) Map view of the sandstone surface colored by height above the “best-fit” plane. (d) Probability density of heights and fitted normal distribution for sandstone surface.

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

Figure 14 Calculated Hurst exponent for different sizes of the measurement area in the natural rock samples. Circles: individual measurements, dashed lines: average of all measurements per sample.

Figure 15 Estimated JRC values calculated based on Eq. (5) for fracture surfaces of the sandstone and limestone specimen. Black: limestone, red: sandstone. Open symbols: profiles taken parallel to shortening direction, filled symbols: profiles perpendicular to shortening direction. Small horizontal offset between data points in each group added for better visibility of individual error bars.

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

5 Discussion

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₂ 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₂ is the rms of the Z₂ 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 σ₂=6 MPa and σ₃=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).

Table 1 Roughness properties for surfaces generated in numerical simulations of triaxial compression tests at σ₂=6 MPa, σ₃=0 using different particle size ranges for the DEM material.

In numerical models, there is a slightly higher anisotropy in the models with transversely isotropic confinement (σ₂=σ₃) of up to 8 % difference in JRC between the directions, whereas in the models with σ₂≠σ₃ the difference is less than 3 % in all cases. In the rock samples, which are also deformed under conditions where σ₂=σ₃, 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 (σ₂≠σ₃).

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 σ₂=7.5 MPa, σ₃=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 (σ₂≠σ₃). 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 σ₂≠σ₃, but that there is no such constraint if σ₂=σ₃. 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.

Figure 16 Comparison of the height–height correlation functions of a numerical fracture surface based on raw data (crosses) and corrected for the influence of the particle-scale roughness (circles). Lines are power-law fits used to calculate the Hurst exponents for raw (continuous line, H≈0.49) and corrected data (dashed line H≈0.53).

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.

Figure 17 Relation between joint roughness coefficient and Hurst exponent of surfaces from numerical models (squares) and the limestone sample (cross). Data points show averages for groups of surfaces generated under the same stress conditions. Error bars on the limestone data shows anisotropy of the JRC. The dashed line shows the relation proposed by Ficker (2017), Eq. (22).

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 σ₃≈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 σ₃=0–15 MPa.

6 Conclusions

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.