Intermittency of Rock Fractured Surfaces: A Power Law

Abstract

Roughness of rock fractured surfaces is one of the most important factors controlling fluid flow in rock masses. Roughness quantification is of prime importance for modelling the flow of ground waters as well as reservoir fluid mechanics. In this study, with the aid of high-resolution 3D X-ray CT scanning and image processing techniques, the roughness of four different rock types is reconstructed with a resolution of 16.5 microns. Moreover, the correlation and structure functions are used to analyse height fluctuations as well as statistical intermittency of the studied rock fractured surfaces. It is observed that at length scales smaller than a critical length scale, fractures surfaces are correlated and show multifractality. Monofractals are neither intermittent nor correlated; hence, a meaningful link between statistical intermittency and the correlation function of multifractals is expected. However, a model that considers this relationship and predicts multifractal spectra of disordered systems is still missing. A simple power law that can exactly forecast the multiscaling spectrum of rock fracture process zone is being introduced. It is explained how the exponent of this power function λ_i is related to the crossover length of correlation function ξ, and how this critical length scale can be objectively identified.

1. Introduction

Quantitative roughness analysis of rock materials is not only the key to understand mysteries associated with nonlinear inelastic fracture mechanics at small enough length scales, but also a powerful and applicable tool in many rock engineering fields [1,2,3].

The main goal of this work is to establish a relationship between the correlation and intermittency of rock fractured surfaces. This link not only broadens our horizons about the statistical physics of multifractal phenomena but also is a promising approach in order to objectively determine the crossover length (or time) $\xi$ of physical phenomena. $\xi$ is a very important parameter in studying the phase transition of natural phenomena [4]. If a fracture is considered as a phase transition [5], $\xi$ is a cross-over length scale that shows a transition from continuous damage percolation at $\delta r\ll \xi$ to the first order at $\delta r\gg \xi$ where material can be considered as linear elastic [6], where $\delta r$ denotes separation or length scale. This phase transition is not very sharp in quasi-brittle rock materials and takes place in a range of length scales because of the mixed first-order and continuous character [7]. Therefore, determining the true $\xi$ as the statistical critical point at which the phase transition takes place, which is reminiscent of the effective length of the fracture process zone $l_{pz}$ [8], is of prime importance and the key to modelling fracture properties of quasi-brittle materials [9,10]. Considering this link, the height–height correlation function [11] or structure function [12] and correlation (auto-covariance) function are integrated in order to model the multiscale spectrum of fracture process zone (FPZ).

2. Materials and Methods

For this work, four different rock types including sandstone, marble and two granites with different grain sizes have been studied. The fractured surfaces are generated by performing notched semi-circular bending tests and under controlled conditions. The topography of fractured surfaces, illustrated in Figure 1, has been reconstructed from 3D X-ray computed tomography data with a spatial resolution $d$ of about $16.5 \mathsf{\mu }\mathrm{m}$. This method is superior to contact methods affected by the tip geometry of probes [13]. The average grain size of the studied rocks is quantified by some 2D slices of tomographic images and ranges from about 0.1 mm for the sandstone to around 1 mm for the marble. It is notable that before any post-processing on fractured surfaces, the global slope in the direction of crack propagation is removed by applying the gradient method such that the average height at the first and last lines of the height map in the propagation direction is the same.

Figure 2 is schematically illustrating the adopted methodology for modelling the multiscaling spectrum of FPZ of studied rocks. First, an initial guess of cross-over length $\xi$ is made by extrapolating the slope of the correlation function ${\lambda }_c$ calculated at some length scales $ϵ$ to $C_{ϵ}\left(\delta r\right)=0$ at which $\delta r=\xi$. Then, the generalized Hurst exponent $H\left(q\right)$ is calculated for both monoaffine ($\delta r>\xi$) and multiaffine ($\delta r<\xi$) regimes. As opposed to monoaffine regime, at multiaffine one $H\left(q\right)$ nontrivially changes with $q$. This intermittency can be modelled with a perfect power law characterised by an exponent ${\lambda }_i$. If ${\lambda }_i=c\ast {\lambda }_c$, it means the initial estimate of $\xi$ was correct and both the true $\xi$ and ${\lambda }_c$ are correctly identified; otherwise, this estimate should be adjusted until ${\lambda }_i=c\ast {\lambda }_c$. $c$ is a constant, and based on experimental observations in this study, its value is considered to be related to the logarithmic mean ${\log }_{e}10/{\log }_{10}10\approx 2.3$. After the first iteration, the point of convergence of multiscale correlation functions and ${\lambda }_i$ are used to modify ${\lambda }_c$, and optimize the initial guess of $\xi$, accordingly.

3. Results and Discussions

3.1. Roughness Correlation

These fractured surfaces exhibit anisotropic scaling properties. In order to identify the effective $\xi$ by means of the correlation function, height variations of fractured surfaces $\Delta h\left(ϵ\right)$ at some length scales $ϵ$ are made statistically isotropic. ${\omega }_{ϵ}\left(X\right)$ the operator has been used for such purpose [14]:

$${\omega }_{ϵ}\left(X\right)=\frac{1}{2}\log \left(⟨\Delta h{\left(ϵ\right)}^2⟩{}_{\mathsf{\Theta }}\right)-{\mathsf{\Omega }}_{ϵ},$$

(1)

where $\log \left(\Delta h\right)$ denotes the natural logarithm, whose base value is $e$. It measures the average height variations around each point, on a height field $h\left(X\right)$, with its neighbours over a circle of radius $ϵ=d\ast p \left(p=1, 2, 3, 4, 5 px\right)$. $ϵ$ is a product of spatial resolution and the number of pixels used for calculating $\Delta h\left(ϵ\right)$; ${\mathsf{\Omega }}_{ϵ}$ is chosen such that the average of ${\omega }_{ϵ}\left(X\right)$ overall $X$ is zero. The pair length scale $ϵ$ and direction $n$ $\left(ϵ,n\right)=\left\{\left(16.5,8\right),\left(33,12\right),\left(49.5,16\right),\left(66,20\right),\left(82.5,28\right)\right\}$ is used to compute $⟨\Delta h{\left(ϵ\right)}^2⟩{}_{\mathsf{\Theta }}=\frac{1}{n}{{\displaystyle \sum }}_{k=0}^{n-1}{\left[h\left(x_i,z_j\right)-h\left(x_i+ϵ\ast \cos \left(2\pi k/n \right),z_j+ϵ\ast \sin \left(2\pi k/n \right)\right)\right]}^2$ where $\mathsf{\Theta }\in \left[0,2\pi \right) rad$. For instance, at a length scale of $ϵ=16.5 \mathsf{\mu }\mathrm{m}$, 8 different data points are used to calculate the height difference for each point. $x_i$ and $z_j$ are the coordinates of a point on $h\left(X\right)$ in crack propagation and crack front directions, respectively. Figure 3 is showing the correlation functions $C_{ϵ}\left(\delta r\right)={⟨{\omega }_{ϵ}\left(X\right){\omega }_{ϵ}\left(X+\delta r\right)⟩}_{\theta }$ of ${\omega }_{ϵ}$ fields averaged over $4$ directions $\theta \in \left[0,\pi \right)=\pi k/4 rad \left(k=0, 1, 2, 3\right)$. The best-fit line of correlation functions $C_{ϵ}\left(\delta r\right)=-{\lambda }_c\log \left(\delta r/\xi \right)+\epsilon$ passing through the point of convergence and the true $\xi$ with the slope ${\lambda }_c={\lambda }_i/c$ for different rock types are shown (under the assumption that $\epsilon =0$). It is notable that ${\lambda }_c$ is dimensionless, since slopes of auto-correlation (dimensionless) and auto-covariance (dimensional) functions are the same. Indeed, $C_{ϵ}\left(\delta r\right)$ indicates scale-dependency of the material’s disorder since ${\omega }_{ϵ}$ quantifies local height variations and removes global slopes caused by macroscopic and dynamic effects.

3.2. Multifractality of Roughness

Following [15], the $q$th root of the $q$th moment of a statistical distribution of the height fluctuations $p\left(\Delta h\right)$, known as $q$th-order structure function $S_q$, is utilized to demonstrate scale-dependency of the roughness of quasi-brittle fractured surfaces, and their transition from monoaffine to multiaffine surfaces at small enough separations:

$$S_q\left(\delta r\right)={{⟨\left|\Delta h\left(\delta r\right)\right|}^q⟩}^{1/q}=⟨{\left|h\left(\mathbf{X}+\delta r\right)-h\left(\mathbf{X}\right)\right|}^q⟩{}_{\mathbf{X}}^{1/q}\propto \delta r^{{\zeta }_q\ast 1/q},$$

(2)

where angular brackets denote the ensemble average.

The concept of multiaffine fractals was introduced in the study of roughness analysis of growing surfaces [16]. FPZ is a multiaffine phenomenon showing different scaling properties in different directions $\theta \in \left[0,\pi \right) rad$ at different separations. However, it shall be shown that its average scaling properties over all directions show a multifractal process. In order to compute the multiscaling spectrum of the fracture process zone, $S_q\left(\theta ,\delta r\right)$ has been calculated over admissible coordinates $\left(x{o}_i,z{o}_j\right)$. $x{o}_i$ and $z{o}_j$ are the coordinates of a point on a height field with zero average, $⟨h\left(X\right)⟩=0$, in crack propagation and crack front directions, respectively. It is a very time-consuming process to compute many moments of structure function over all admissible coordinates and directions with small intervals for large data sets in order of million data points. $m=12$ moments, $q=\left\{0.2, 0.4, 0.6, 0.8, 1, 1.5, 2, 2.5, 3, 4, 5, 6\right\}$, $n=4$ directions, and $\theta =\pi k/n rad \left(k=0,1,\dots ,n-1\right)$ have been selected, which are accurate enough for analysing the intermittency of the studied rocks. Moreover, an interval $l=2 px$ for calculating $S_q\left(\theta ,\delta r\right)$ is used. Thus, $\Delta h\left(\delta r\right)$ has been calculated between the following admissible coordinates $h\left(x{o}_{il},z{o}_{jl}+\delta r\right)-h\left(x{o}_{il},z{o}_{jl}\right)$ for selected moment and angle sets. It is notable that for a few smaller subsets $S_q\left(\theta ,\delta r\right)$ have been calculated among both selected and all admissible coordinates $h\left(x{o}_i,z{o}_j+\delta r\right)-h\left(x{o}_i,z{o}_j\right)$ in different directions $n=\left\{4, 12, 36\right\}$, and differences between calculated average scaling properties $⟨S_q{\left(\delta r\right)}_{\theta }⟩$ were negligible. Admissible coordinates would be different considering separation and direction. Maximum separation has always been about half of the minimum dimension of the studied surfaces to compute moments at different separations over enough and similar data points. Overall, 34 to 37 different separations $\delta r\in \left[16.5,13200\right] \mathsf{\mu }\mathrm{m}$ have been selected considering the dimensions of studied surfaces. Multiscaling spectra of FPZ of studied rocks are presented in Figure 4.

3.3. Determining the Cut-Off Length

It is observed that heart rate time series display different multiscaling exponents because of long-range correlations [17]. The BDM model is proposed to provide a link between multifractality and long-range correlation of financial time series [18]:

$$H\left(q\right)=H-\left(q-1\right)\frac{{\lambda }_c}{2},$$

(3)

in which $H\equiv {\zeta }_1$, and ${\lambda }_c$ is the slope of the correlation function of a multifractal process. This means that $H\left(q\right)=H$ if there is no correlation, i.e., ${\lambda }_c=0$, which is the case for monofractals. Otherwise, $H\left(q\right)$ decrease linearly with $q$ for ${\lambda }_c>0$. The BDM model has been used to analyse the multifractal behaviour of FPZ [14]. However, this linear model could not provide a good estimation of such non-linear multifractal spectra of the studied rocks. In this study, an experimental power law is introduced that can exactly predict the statistically isotropic multiscaling spectrum of rock FPZ:

$$H\left(q\right)=H\left(q^{-{\lambda }_i}\right),$$

(4)

in which ${\lambda }_i$ is the intermittency (slope of the multiscaling spectrum) of a multifractal process. This model is derived from experimental observation of the intermittent behaviour of the multiscaling spectrum of FPZ of the studied rock materials.

By applying a log–log transformation, Equation (4) can be read as $\log H\left(q\right)=-{\lambda }_i\log q+\log H$. From Figure 4, it is indisputable that this log–log transformation is showing a perfect linear relationship that is unique to power functions. ${\lambda }_i$ computed from this relationship has been used to adjust ${\lambda }_c$ and optimize the initial guess of $\xi$. From Figure 3, it can be clearly seen that apparent correlation functions computed at different scales meet each other at a particular point, which is the true correlation function and reflects the intrinsic length of the material disorder. The coordinates of this point are the key to determine the true $\xi$ by modifying ${\lambda }_c$ based on computing ${\lambda }_i$. By repeating the loop shown in Figure 1 a few times (from 5 to 8 iterations in this study) both the intermittency of the multifractal process and the true cross-over length of the correlation function can be objectively identified. A comparison among experimental data, the BDM model and the proposed model is made in Figure 5, which shows the success of the proposed model in predicting multiscaling spectra of the rock fractured surfaces at small enough length scales.

If ${\lambda }_i=c\ast {\lambda }_c$, then $1-\frac{\log \left({\zeta }_q/H\right)}{\log \left(q\right)}=-c\ast \frac{C_{ϵ}^{*}\left(\delta r\right)-\epsilon }{\log \left(\delta r/\xi \right)}$, and this equation can be rearranged as follows:

$$\frac{\log \left({\zeta }_q/H\right)}{\log \left(q\right)}-\frac{c\ast C_{ϵ}^{*}\left(\delta r\right)}{\log \left(\delta r/\xi \right)}=1-\frac{c\ast \epsilon }{\log \left(\delta r/\xi \right)},$$

(5)

where $C_{ϵ}^{*}\left(\delta r\right)$ is an auto-correlation function that is normalized by the variance of ${\omega }_{ϵ}$ field. In the monoaffine regime, ${\zeta }_q/H=q$, there is no correlation $C_{ϵ}^{*}\left(\delta r\right)\approx 0$, and Equation (5) is read as $1-0=1$. In the transition zone, the correlation from $\epsilon$ at $\delta r=\xi$ would tend to zero as $\delta r$ tends to infinity. In the multiaffine regime, ${\zeta }_q/H=q^{1-{\lambda }_i}$, there is a correlation of $0<C_{ϵ}^{*}\left(\delta r\right)-\epsilon <1$, and Equation (5) read as $1-{\lambda }_i+c\ast {\lambda }_c=1$.

3.4. Limitations and Future Work

The question that remains open here is as follows: what is the relationship between the slope of the correlation function and the intermittency? To answer this question, finding a meaningful link between the correlation function and the structure function is the first step. Two kinds of these links that can trigger some interesting future research are provided as follows.

The power spectrum, on the one hand, is the Fourier transform of the auto-covariance function of a wide-sense stationary random process (Wiener–Khinchin theorem). Height variations in the monofractal regime can be considered as a fractional Brownian motion with stationary increments and their spectral density $S\left(f\right)$ scales with the frequency $f$ as a power law: $S\left(f\right)\propto f^{-\beta }$ where $\beta =2H\left(2\right)+1$ [19]. On the other hand, the spectral density is proportional to the second moment of distribution of height variations $⟨{\left|\Delta h\left(\delta r\right)\right|}^2⟩\propto \delta r^{\alpha }$ where $\alpha =2H\left(2\right)=\beta -1$. For monofractals a relationship between $H\left(q\right)=H$ and $\beta$ can be defined as $D_q=q+1-H=q+\frac{3-\beta }{2}$. However, when it comes to multifractals with nonstationary increments, this link between the correlation function and structure function needs some modifications. In other words, like the generalized fractal dimension and Hurst exponent, a generalization for spectral analysis is required. Wigner–Ville spectral analysis is introduced as the unique generalized spectrum for spectral analysis of non-stationary processes [20]. There is also a statistical relationship between the structure function and correlation function: $⟨{\left|\Delta h\left(\delta r\right)\right|}^2⟩=2\left({\sigma }^2-C_{ϵ}\left(\delta r\right)\right)$, where ${\sigma }^2=⟨h{\left(X\right)}^2⟩\approx ⟨h{\left(X+\delta r\right)}^2⟩$ is the variance of a height map. This equation is working very well for the second moment, but the question is how it can be expanded for other moments. Moreover, this relationship ignores the effect of local averaging on the height map (to compute ${\omega }_{ϵ}$ field). Along with this theoretical issue, it seems the coefficient $c$ plays a very important role in understanding multifractal phenomena and can be used to estimate the crossover length of multiphase phenomena. Further measurements in different fields are required to have a better comprehension of the proposed model in this study. The outcomes of this research can pave the way to deciphering multiscaling features of quasi-brittle fractured surfaces and model the roughness of the fractured rocks.

4. Conclusions

In this study, the intermittency of the roughness of fractured surfaces of different rock types is quantified. Based on the studied fractured surfaces the intermittency can be modelled by perfect power laws that can be further used to analyse multiscale properties of the roughness of fractured surfaces at length scales smaller than the cut-off length, an important topic in fluid mechanics. Moreover, a repetitive loop is introduced to determine the cut-off length, which is the length scale at which a transition from multi- to monofractality of the fracture surfaces can be observed, by making a connection between the slope of intermittency and slope of the correlation function of the same surface. Finally, some directions for future research are introduced that can further enhance our understanding of the multiscaling features of rough surfaces.