A Fractal-Enhanced Mohr–Coulomb Model for Strength Prediction in Rough Rock Discontinuities

1. Introduction

In many engineering disciplines dealing with natural materials, the assessment of geotechnical stability is fundamental for many design decisions. Given the predisposition of structures such as rock slopes, tunnels, and foundations, the key failure mechanism can be considered as material damage, often represented as a generalized failure boundary in the stress space. Several material failure criteria have been formulated over the past century. The Mohr–Coulomb criterion is one of the most fundamental formulations underpinning practical engineering applications owing to its clear physical interpretation and straightforward mathematical representation [1,2,3,4,5,6,7]. The shear resistance can be viewed as a linear combination of two components, an intrinsic cohesive part and a frictional part, provided that the material properties are uniform and the failure surface can be planar. This theory can provide good results, but it can be imprecise because of the complex nature of the constitutive responses of natural geomaterials, which are commonly nonhomogeneous and discontinuous[8,9,10]. Real rock joints, faults, and stratification planes have complex surface roughness and waviness that characterize their morphology and bearing capacity. This feature is not captured by the classical Mohr–Coulomb theory because surface irregularities are a function of physical reality rather than material properties [11,12,13,14]. Although roughness can be characterized in practice by some means, these treatments are chiefly empirical and are based on indirect visual determinations of profile measurements.

Fractal dimension represents the single-parameter scaling law of fractal geometry, which has been extensively studied for natural surfaces of harsh scales. Unlike Euclidean geometry, fractal dimension reflects the self-affine behavior of rock fractures and invariant features at different magnifications and sampling resolutions [15,16,17,18,19,20]. The fractal dimension has been rigorously validated in previous studies, with a strong statistical association with experimental shear strength [21,22,23,24,25]. Although fractal dimension has always been used as an inverse parameter in regression models and classifying algorithms in existing research, it has not been extensively used as a direct mechanical input parameter embedded within the functional form of a failure criterion. This implies that mechanical laws can be fundamentally refined and extended in new ways. In this context, having a real, quantifiable, and mechanical geometric space represents a remarkable opportunity for geomechanical modeling. To ensure that the fractal dimension serves as a reliable and objective mechanical input parameter, a robust measurement method was applied. Topographic data contribute to box counting and power spectral density analyses of well-known procedures. The two methods were used independently to provide cross-validation of the results; the first was based on spatial coverage, and the second on frequency-space characteristics. The contradictions between them were minimized via appropriate pre-processing and selection of the scale range, thereby increasing the reliability of the group parameter estimation. Through this dual method validation, there was higher certainty in the volumetric results and a lower propagation risk of bias and errors. The accuracy level provided for the fractal measurements is guaranteed to allow the scaling of such parameters into the physical harm criterion.

This study aims to develop a fractal-enhanced Mohr–Coulomb failure criterion by directly embedding the surface fractal dimension DD as a constitutive parameter that governs the effective cohesion $c\left(D_S\right)$ and internal friction angle $\phi \left(D_S\right)$. By transforming $D_S$ from a descriptive geometric statistic into an active mechanical variable, this study establishes a direct predictive link between surface topography and macroscopic strength behavior.

2. Methodology

This study employs a systematic approach to develop an enhanced Mohr-Coulomb failure criterion through direct integration of the surface fractal dimension[26]. The resulting fractal dimension values were then used to formulate scale-dependent expressions for the cohesion ($c\text{'}$) and friction angle ($\phi \text{'}$), which formed the basis of the proposed strength model.

2.1. Mathematical Description of Fractal Dimension ($D_{\mathrm{S}}$)

A box-counting method[27] was included in fractal geometry to calculate the fractal dimension. The idea is to cover an object of irregular shape with smaller boxes (or grids) and count the number $\mathsf{{\rm N}}\left(\epsilon \right)$ of boxes required to cover the object as a function of the box size $\epsilon$.

The relation can be described as:

$$\mathsf{{\rm N}}\left(\epsilon \right)=C\cdot {\epsilon }^{{-\mathrm{D}}_{\mathrm{S}}}$$

(1)

where:

• $\mathsf{{\rm N}}\left(\epsilon \right)$ is the number of boxes,

• $C$ is a constant depending on the shape,

• $D_{\mathrm{S}}$ is the fractal dimension of the object.

Taking the natural logarithm of both sides yields a linear relationship:

$$\ln \mathsf{{\rm N}}\left(\epsilon \right)=\ln C-D_{\mathrm{S}}\ln \epsilon$$

(2)

This equation mirrors a straight line, where the slope corresponds to the fractal dimension ($D_{\mathrm{S}}$)

• The expression can be further simplified into a standard linear form: .

  $$y=mx+b$$

  (3)
  where:
  • y = ln N(ε).
  • x = ln (ε),
  • m = −D_s (negative slope).
  • b = ln C (the y-intercept).
  The fractal dimension $D_{\mathrm{S}}$was then computed as the negative slope of the line fitted to the experimental data points plotted in a log-log space.

2.2. Power Spectral Density Method: Frequency Domain Analysis

The power spectral density (PSD)[28,29,30] analyzes surface topography based on its spatial frequency content. Given a surface height function, $h(x)$, its Fourier transform yields PSD $\mathrm{{\rm P}}\left(k\right)$, where $k$ is the spatial frequency (wavenumber). For a fractal surface, the PSD follows a power-law decay in the frequency domain:

$$\mathrm{{\rm P}}\left(k\right)=\mathrm{A}\cdot k^{-\left(5-2D_S\right)}$$

(4)

where $\mathrm{A}$ is a scaling constant related to the surface amplitude. Taking the logarithm of both sides converts Equation (4) into a linear form, as follows:

$$\ln \mathrm{{\rm P}}\left(k\right)=\ln \mathrm{A}-\left(5-2D_{\mathrm{S}}\right)\ln k$$

(5)

The slope$m$of the $\ln \mathrm{{\rm P}}\left(k\right)$versus $\ln k$ plot directly yields the fractal dimension:

$$m=-\left(5-2D_{\mathrm{S}}\right)$$

(6)

$$\begin{array}{l}{D}_{\mathrm{S}}={\displaystyle \frac{5+m}{2}}\end{array}$$

(7)

The PSD method offers superior noise resistance compared with box counting, as frequency-domain analysis naturally filters high-frequency measurement artifacts. Both methods were applied in this study using cross-validation to ensure robustness.

2.3. Classical Mohr-Coulomb theory

The shear strength ${\tau }_f$at failure is defined as:

$${\tau }_f=c\text{'}+\sigma \text{'}\tan \phi \text{'}$$

(8)

• c^′ is the effective cohesion,
• σ^′ is the effective normal stress on the failure plane, and
• φ^′ where denotes the effective friction angle. Under triaxial testing conditions, the failure criterion was expressed in terms of the major and minor principal stresses. We begin with the geometric relationship of stress transformation. The radius of the Mohr circle at failure is

where:

$$R={\displaystyle \frac{{\sigma }_1^{\text{'}}-{\sigma }_3^{\text{'}}}{2}}$$

(9)

The center of the circle is located at

$$C={\displaystyle \frac{{\sigma }_1^{\text{'}}+{\sigma }_3^{\text{'}}}{2}}$$

(10)

The failure plane makes an angle ${\theta }_f$ with the major principal plane, where

$${\theta }_f={45}^o+{\displaystyle \frac{\phi \text{'}}{2}}$$

(11)

At failure, the Mohr circle is tangent to the failure envelope. The geometric condition requires the following:

$$R\sin \phi \text{'}=C\sin \phi \text{'}+c\text{'}\cos \phi \text{'}$$

(12)

Substituting the expressions for $R$ and $C$:

$${\displaystyle \frac{{\sigma }_1^{\text{'}}-{\sigma }_3^{\text{'}}}{2}}\sin \phi \text{'}={\displaystyle \frac{{\sigma }_1^{\text{'}}+{\sigma }_3^{\text{'}}}{2}}\sin \phi \text{'}+c\text{'}\cos \phi \text{'}$$

(13)

Multiplying both sides by 2 and expanding:

$$\left({\sigma }_1^{\text{'}}-{\sigma }_3^{\text{'}}\right)\sin \phi \text{'}=\left({\sigma }_1^{\text{'}}+{\sigma }_3^{\text{'}}\right)\sin \phi \text{'}+2c\text{'}\cos \phi \text{'}$$

(14)

Collecting terms with ${\sigma }_1^{\text{'}}$ and ${\sigma }_3^{\text{'}}$ :

$${\sigma }_1^{\text{'}}\sin \phi \text{'}-{\sigma }_3^{\text{'}}\sin \phi \text{'}={\sigma }_1^{\text{'}}\sin \phi \text{'}+{\sigma }_3^{\text{'}}\sin \phi \text{'}+2c\text{'}\cos \phi \text{'}$$

(15)

$$-2{\sigma }_3^{\text{'}}\sin \phi \text{'}=2c\text{'}\cos \phi \text{'}$$

(16)

$${\sigma }_1^{\text{'}}-{\sigma }_1^{\text{'}}\sin \phi \text{'}={\sigma }_3^{\text{'}}+{\sigma }_3^{\text{'}}\sin \phi \text{'}+2c\text{'}\cos \phi \text{'}$$

(17)

Factoring and dividing (17) by $1-\sin \phi \text{'}$, we have

$${\sigma }_1^{\text{'}}={\sigma }_3^{\text{'}}{\displaystyle \frac{1+\sin \phi \text{'}}{1-\sin \phi \text{'}}}+2c\text{'}{\displaystyle \frac{\cos \phi \text{'}}{1-\sin \phi \text{'}}}$$

(18)

Compact the notation with ${{\rm N}}_{\phi }$ , defining the bearing capacity factor:

$${{\rm N}}_{\phi }={\tan }^2\left({\displaystyle \frac{\phi \text{'}}{2}}+{45}^o\right)={\displaystyle \frac{1+\sin \phi \text{'}}{1-\sin \phi \text{'}}}$$

(19)

The classical Mohr-Coulomb criterion becomes:

$${\sigma }_1^{\text{'}}={\sigma }_3^{\text{'}}{{\rm N}}_{\phi }+2c\text{'}\sqrt{{{\rm N}}_{\phi }}$$

(20)

The term ${{\rm N}}_{\phi }$represents the amplification of the confining stress due to friction, whereas $2c\text{'}\sqrt{{{\rm N}}_{\phi }}$ represents the contribution of the cohesive strength. Both terms are assumed to be constant in classical theory. We relax this assumption through fractal enhancement.

2.4. Fractal Enhancement of Mohr-Coulomb Theory

Experimental evidence demonstrates that surface roughness significantly influences mechanical behavior. Rougher surfaces ($D_{\mathrm{S}}$low) exhibit enhanced interlocking, leading to higher peak strength and greater dilation. Conversely, smoother surfaces ($D_{\mathrm{S}}$ high) result from weathering or particle crushing, producing reduced frictional resistance and easier strain softening. We propose that the classical strength parameters are not material constants, but rather functions of the measurable fractal dimension:

$$c\text{'}=c\text{'}\left(D_{\mathrm{S}}\right)$$

(21)

$$\phi \text{'}=\phi \text{'}\left(D_{\mathrm{S}}\right)$$

(22)

Based on experimental observation across multiple geological classes, the effective internal friction angle $\phi \text{'}$ exhibits approximately linear dependence on fractal dimension:

$$\phi \text{'}\left(D_{\mathrm{S}}\right)={\phi }_0+k_{\phi }\left(D_{\mathrm{S}}-D_0\right)$$

(23)

where:

• ${\phi }_0$ : Reference friction angle at performable fractal dimension $D_0$
• $k_{\phi }$ : Sensitivity coefficient (units: degrees per unit $D_{\mathrm{S}}$)
• $D_0$ : Reference fractal dimension, conventionally set at $D_0=2.0$ for perfectly smooth Euclidean surfaces;

• $D_{\mathrm{S}}$: The fractal dimensions of an actual material surface were measured using box-counting or PSD methods.

Effective cohesion exhibits a more complex non-linear dependence on surface roughness owing to the competing effects of aspect interlocking and contact-area evolution.

$$c\text{'}\left(D_{\mathrm{S}}\right)=c_0+c_1^{\left[k_c\left(D_{\mathrm{S}}-D_0\right)\right]}$$

(24)

where:

• c₀ = Intrinsic cohesion (kPa or MPa)
• c₁ = roughness amplification factor.
• k_c = exponential decay rate.

Substituting the fractal-dependent strength parameters into the classical Mohr-Coulomb criterion (8) yields the following enhanced model:

$${\tau }_f\left(D_{\mathrm{S}}\right)=\left\{c_0+c_1^{\left[k_c\left(D_{\mathrm{S}}-D_0\right)\right]}\right\}+{\sigma }_n^{\text{'}}\tan \left[{\phi }_0+k_{\phi }\left(D_{\mathrm{S}}-D_0\right)\right]$$

(25)

In principal stress space, the modified criterion becomes:

$${\sigma }_1^{\text{'}}={\sigma }_3^{\text{'}}{{\rm N}}_{\phi }\left(D_{\mathrm{S}}\right)+2c\text{'}\left(D_{\mathrm{S}}\right)\sqrt{{{\rm N}}_{\phi }\left(D_{\mathrm{S}}\right)}$$

(26)

where the bearing capacity factor is fracture dependent.

$${{\rm N}}_{\phi }\left(D_{\mathrm{S}}\right)={\tan }^2\left({45}^o+{\displaystyle \frac{\phi \text{'}\left(D_{\mathrm{S}}\right)}{2}}\right)$$

(27)

This formulation naturally captures the evolution of the field surface with microstructural roughness, enabling predictions of strength variation without empirical recalibration for each material type.

2.5. Fractional Calculus Framework

Classical integer-order differential equations assume an instantaneous response and local behavior. However, geomaterials exhibit memory effects, nonlocal interactions, and power-law creep phenomena, which are better described by fractional calculus. The Caputo fractional derivative of order$\alpha \in \left(0,1\right)$ is defined as:

$${}^{C}D_t^{\alpha }f\left(t\right)={\displaystyle \frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)}}{\displaystyle {\int }_0^{\text{'}}\frac{\left(\tau \right)}{{\left(t-\tau \right)}^{\alpha }}}d\tau$$

(28)

where $\mathrm{\Gamma }$ is the gamma function. This operator was reduced to classical derivatives during $\alpha =1$. The Atangana-Baleanu (AB) derivative in the Caputo sense [31,32]overcomes the singularity limitation of classical fractional operators.

$${}^{\mathit{ABC}}D_t^{\alpha }f\left(t\right)={\displaystyle \frac{\mathrm{{\rm B}}\left(\alpha \right)}{1-\alpha }}{\displaystyle {\int }_0^{t}f\text{'}\left(\tau \right)}{\mathrm{{\rm E}}}_{\alpha }\left[-{\displaystyle \frac{\alpha }{1-\alpha }}{\left(t-\tau \right)}^{\alpha }\right]d\tau$$

(29)

where $\mathrm{{\rm B}}\left(\alpha \right)$is the normalization function satisfying:

$\mathrm{{\rm B}}\left(0\right)=\mathrm{{\rm B}}\left(1\right)=1$, which ensures consistency with integer-order limits.

The function ${\mathrm{{\rm E}}}_{\alpha }\left(z\right)$ is given by:

$${\mathrm{{\rm E}}}_{\alpha }\left(z\right)={\displaystyle \sum _{k=0}^{\infty }\frac{z^k}{\mathrm{\Gamma }\left(\alpha k+1\right)}}$$

(30)

The Laplace transform of the Caputo derivative is a powerful analytical tool:

$$\wp \left\{{}^{C}D_t^{\alpha }f\left(t\right)\right\}=s^{\alpha }F\left(s\right)-s^{\alpha -1}f\left(0\right)$$

(31)

For the AB derivative, the Laplace transform is:

$$\wp \left\{{}^{\mathit{ABC}}D_t^{\alpha }f\left(t\right)\right\}={\displaystyle \frac{B\left(\alpha \right)}{1-\alpha }}\cdot {\displaystyle \frac{s^{\alpha }F\left(s\right)-s^{\alpha -1}f\left(0\right)}{s^{\alpha }+\frac{\alpha }{1-\alpha }}}$$

(32)

The fractal-enhanced Mohr-Coulomb model was formulated as a fractional-order stress evolution equation coupling fractal geometry with material memory, as follows:

$${}^{\mathit{ABC}}D_t^{\alpha }\sigma \left(t\right)+\beta \sigma \left(t\right)=\gamma c\text{'}\left(D_{\mathrm{S}}\right)+\delta \phi \text{'}\left(D_{\mathrm{S}}\right)\epsilon \left(t\right)$$

(33)

where:

• σ(t) = Time-Dependent Stress,
• α = fractional order 0<α≤1,
• β = damping coefficient,
• γ, δ = coupling constants,
• ε(t) = strain history.

By applying the Laplace transform to the governing equation with zero initial conditions, we obtain the following general solution:

$${\displaystyle \frac{B\left(\alpha \right)}{1-\alpha }}\cdot {\displaystyle \frac{s^{\alpha }{\displaystyle \sum \left(s\right)}}{s^{\alpha }+\frac{\alpha }{1-\alpha }}}+\beta {\displaystyle \sum \left(s\right)}=\gamma c\text{'}\left(D_{\mathrm{S}}\right){\displaystyle \frac{1}{s}}+\delta \phi \text{'}\left(D_{\mathrm{S}}\right)\mathrm{{\rm E}}\left(s\right)$$

(34)

: ${\displaystyle \sum \left(s\right)}=\wp \left\{\sigma \left(t\right)\right\}$, $\mathrm{{\rm E}}\left(s\right)=\wp \left\{\epsilon \left(t\right)\right\}$. Multiplying by $\left(1-\alpha \right)$ and rearranging the denominator terms:

$${\displaystyle \sum \left(s\right)}\left[\mathrm{{\rm B}}\left(\alpha \right)s^{\alpha }+\beta \left(1-\alpha \right)\left(s^{\alpha }+{\displaystyle \frac{\alpha }{1-\alpha }}\right)\right]=\left[\gamma c\text{'}\left(D_{\mathrm{S}}\right){\displaystyle \frac{1}{s}}+\delta \phi \text{'}\left(D_{\mathrm{S}}\right)\mathrm{{\rm E}}\left(S\right)\right]\left(1-\alpha \right)\left(S^{\alpha }+{\displaystyle \frac{\alpha }{1-\alpha }}\right)$$

(35)

Collect terms and solve for ${\displaystyle \sum \left(s\right)}$:

$${\displaystyle \sum \left(s\right)}={\displaystyle \frac{\left(1-\alpha \right)\left(s^{\alpha }+\frac{\alpha }{1-\alpha }\right)}{\mathrm{{\rm B}}\left(\alpha \right)s^{\alpha }+\beta \left(1-\alpha \right)s^{\alpha }+\beta \alpha }}\left[{\displaystyle \frac{\gamma c\text{'}\left(D_{\mathrm{S}}\right)}{s}}+\delta \phi \text{'}\left(D_{\mathrm{S}}\right)\mathrm{{\rm E}}\left(s\right)\right]$$

(36)

Defining the characteristic polynomial:

$$\mathrm{{\rm P}}\left(s\right)=\left[\mathrm{{\rm B}}\left(\alpha \right)+\beta \left(1-\alpha \right)\right]s^{\alpha }+\beta \alpha$$

(37)

The general solution in the Laplace domain becomes:

$${\displaystyle \sum \left(s\right)}={\displaystyle \frac{\left(1-\alpha \right)s^{\alpha }+\alpha }{\mathrm{{\rm P}}\left(s\right)}}\left[{\displaystyle \frac{\gamma c\text{'}\left(D_{\mathrm{S}}\right)}{s}}+\delta \phi \text{'}\left(D_{\mathrm{S}}\right)\mathrm{{\rm E}}\left(s\right)\right]$$

(38)

In order to obtain the general solution in the time domain, we apply the inverse Laplace transform using the convolution theorem and Mittag-Leffler function[33] properties:

$$\sigma \left(t\right)=\gamma c\text{'}\left(D_{\mathrm{S}}\right){\displaystyle {\int }_0^t\mathrm{{\rm K}}\left(t-\tau \right)}d\tau +\delta \phi \text{'}\left(D_{\mathrm{S}}\right){\displaystyle {\int }_0^t\mathrm{{\rm K}}\left(t-\tau \right)}\epsilon \left(\tau \right)d\tau$$

(39)

Where the kernel function $\mathrm{{\rm K}}\left(t\right)$is given by:

$$\mathrm{{\rm K}}\left(t\right)=t^{\alpha -1}{\mathrm{{\rm E}}}_{\alpha ,\alpha }\left({\displaystyle \frac{{\beta }_{\alpha }}{\beta \left(\alpha \right)+\beta \left(1-\alpha \right)}}{t}^{\alpha }\right)$$

(40)

And

$${\mathrm{{\rm E}}}_{\alpha ,\beta }\left(z\right)={\displaystyle \sum _{k=0}^{\infty }\frac{z^k}{\mathrm{\Gamma }\left(\alpha k+\beta \right)}}$$

(41)

2.6. Existence and Uniqueness via Banach Fixed-Point Theorem

The existence and uniqueness of the solution to the fractional constitutive equation is established (41) using the Banach fixed-point theorem[34]. Consider the integral form derived from the Caputo fractional derivative:

$$\sigma \left(t\right)={\sigma }_0+{\displaystyle {\int }_0^t\mathrm{{\rm K}}\left(t-z\right)}F\left(\sigma \left(\tau \right),D_S\right)d\tau$$

(42)

where:

• σ₀ is the initial stress,
• $\mathrm{{\rm K}}\left(t\right)={\displaystyle \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}}{t}^{\alpha -1}$is the memory kernel,
• $F\left(\sigma ,D_S\right)=-\beta \sigma +\gamma c\text{'}\left(D_{\mathrm{S}}\right)+\delta \phi \text{'}\left(D_{\mathrm{S}}\right)\epsilon \left(t\right)$ is a nonlinear forcing term that is assumed to be continuous in $t$, and Lipschitz continuous in $\sigma$.

The memory kernel is given explicitly by:

$$\mathrm{{\rm K}}\left(t\right)={\displaystyle \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}}{t}^{\alpha -1}$$

(43)

Let ${\rm X}=C\left[0,{\rm T}\right]$ denote the Banach space of continuous functions on $\left[0,{\rm T}\right]$equipped with the supremum norm:

$${‖\sigma ‖}_{\infty }=\underset{t\in \left[0,{\rm T}\right]}{\max }\left|\sigma \left(t\right)\right|$$

(44)

Define the solution operator ${\rm T}:{\rm X}\to {\rm X}$ by:

$$\left({\rm T}\sigma \right)\left(t\right)={\sigma }_0+{\displaystyle {\int }_0^t\mathrm{{\rm K}}\left(t-\tau \right)F\left(\sigma \left(\tau \right),D_S\right)}d\tau$$

Assume that $F\left(\sigma ,D_S\right)$ is Lipschitz continuous in $\sigma$ with constant $\beta >0$ , i.e.,

$$\left|F\left({\sigma }_1,D_S\right)-F\left({\sigma }_2,D_S\right)\right|\le \beta \left|{\sigma }_1-{\sigma }_2\right|,\forall {\sigma }_1,{\sigma }_2\in ℝ$$

. For any ${\sigma }_1,{\sigma }_2\in {\rm X}$ , we have:

$$\left|\left({\rm T}{\sigma }_1\right)\left(t\right)-\left({\rm T}{\sigma }_2\right)\left(t\right)\right|\le {\displaystyle {\int }_0^t\left|\mathrm{{\rm K}}\left(t-\tau \right)\right|\cdot \left|F\left({\sigma }_1\left(\tau \right),D_S\right)\right|}d\tau \le \beta {\displaystyle {\int }_0^t\left|\mathrm{{\rm K}}\left(t-\tau \right)\right|\cdot \left|\left({\sigma }_1\left(\tau \right)-{\sigma }_2\left(\tau \right)\right)\right|}d\tau$$

Taking the supremum over t ∈ [0,T], we obtain:

$${‖{\rm T}{\sigma }_1-{\rm T}{\sigma }_2‖}_{\infty }\le \beta \left(\underset{t\in \left[0,{\rm T}\right]}{\sup }{\displaystyle {\int }_0^t\left|\mathrm{{\rm K}}\left(t-\tau \right)\right|}\right)d\tau ‖{\sigma }_1-{\sigma }_2‖$$

(45)

Since $\mathrm{{\rm K}}\left(t\right)={\displaystyle \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}}{t}^{\alpha -1}$, we compute:

$${\displaystyle {\int }_0^t\left|\mathrm{{\rm K}}\left(t-\tau \right)\right|}d\tau ={\displaystyle {\int }_0^t\frac{1}{\mathrm{\Gamma }\left(\alpha \right)}}{\left(t-\tau \right)}^{\alpha -1}d\tau ={\displaystyle \frac{1}{\mathrm{\Gamma }\left(\alpha \right)}}{\displaystyle {\int }_0^t{S}^{\alpha -1}}dS={\displaystyle \frac{t^{\alpha }}{\mathrm{\Gamma }\left(\alpha \right)\alpha }}={\displaystyle \frac{t^{\alpha }}{\mathrm{\Gamma }\left(\alpha +1\right)}}$$

Therefore:

$$\underset{t\in \left[0,{\rm T}\right]}{\sup }{\displaystyle {\int }_0^t\left|\mathrm{{\rm K}}\left(t-\tau \right)\right|}d\tau ={\displaystyle \frac{{{\rm T}}^{\alpha }}{\mathrm{\Gamma }\left(\alpha +1\right)}}$$

and so:

$${‖{\rm T}{\sigma }_1-{\rm T}{\sigma }_2‖}_{\infty }\le \left(\beta \cdot {\displaystyle \frac{{{\rm T}}^{\alpha }}{\left(\alpha +1\right)}}\right){‖{\sigma }_1-{\sigma }_2‖}_{\infty }$$

(46)

Define:

$q=\beta \cdot {\displaystyle \frac{{{\rm T}}^{\alpha }}{\left(\alpha +1\right)}}$$q<1$ , ${\rm T}$ is a contraction mapping of ${\rm X}$. This condition holds for sufficiently small ${\rm T}>0$ , specifically:

$${\rm T}<\left[{\displaystyle \frac{\mathrm{\Gamma }\left(\alpha +1\right)}{\beta }}\right]1/\alpha$$

Theorem (existence and uniqueness): Under the above assumptions, operator ${\rm T}$ admits a unique fixed-point ${\sigma }^{\ast }\in C\left[0,{\rm T}\right]$such that ${\rm T}{\sigma }^{\ast }={\sigma }^{\ast }$. This fixed point is the unique continuous solution to the integral equation (42) and hence to the fractional constitutive law (41) interpreted in the Caputo sense.

Proof sketch: Using the Banach fixed-point theorem, the iterative sequence $\left\{{\sigma }_n\right\}$ is defined as:

$${\sigma }_{n+1}\left(t\right)={\sigma }_0+{\displaystyle {\int }_0^t\mathrm{{\rm K}}\left(t-\tau \right)F\left({\sigma }_n\left(\tau \right),D_S\right)}d\tau$$

converges uniformly on $\left[0,{\rm T}\right]$ to the unique solution ${\sigma }^{\ast }\left(t\right)$. This completes the proof.

2.7. Lyapunov Stability Analysis

Define the Lyapunov functional:

$$V\left(t\right)={\displaystyle \frac{1}{2}}{\sigma }^2\left(t\right)$$

(47)

Computation of the fractional derivative

$${}^{C}D_t^{\alpha }V\left(t\right)=\sigma \left(t\right)\cdot {}^{C}D_t^{\alpha }\sigma \left(t\right)$$

(48)

Substituting the governing equation:

$${}^{C}D_t^{\alpha }V\left(t\right)=\sigma \left(t\right)\left[-\beta \sigma \left(t\right)+\gamma c\text{'}\left(D_{\mathrm{S}}\right)+\delta \phi \text{'}\left(D_{\mathrm{S}}\right)\epsilon \left(t\right)\right]$$

(49)

For bounded loading and positive damping ($\beta >0$), we obtain:

$${}^{C}D_t^{\alpha }V\left(t\right)\le -\beta {\sigma }^2\left(t\right)+\sigma \left(t\right)\cdot C_{bounded}$$

(50)

where: $C_{bounded}=\max \left\{\gamma c\text{'}\left(D_{\mathrm{S}}\right)+\delta \phi \text{'}\left(D_{\mathrm{S}}\right)\epsilon \left(t\right)\right\}$. This implies an exponential-like decay in the fractional sense that ensures stability.

Stability theorem: The solution $\sigma \left(t\right)$is Mittag-Leffler stable. There exists a constant $\lambda >0$and $b>0$ such that:

$$\left|\sigma \left(t\right)\le b\mathrm{{\rm E}}\alpha \left(-\lambda t\alpha \right)\right|$$

(51)

This fractional stability generalizes the exponential stability and is appropriate for systems with memory and nonlocal effects.

3. Results

This paper presents a comprehensive framework for modifying the classical Mohr–Coulomb (MC) failure criterion by incorporating the surface fractal dimension $D_S$ as a state-dependent parameter. The resulting model integrated fractal geometry, fractional calculus, and nonlinear stability theory to predict the strength, deformation, and failure of rock joints under monotonic and cyclic loading. The results demonstrate that the model is not only phenomenologically accurate, but also mathematically well-posed and numerically robust.

3.1. Fractal Dimension Quantification Using Box-Counting and PSD Methods

The fractal dimensions of the rock joint surfaces were quantified using two independent methods: the box-counting method (BCM) and Power Spectral Density (PSD) analysis. The log-log plots in Figure 1a (PSD) and Figure 1b (BCM) exhibit strong linear trends with correlation coefficients exceeding $R^2=0.96$ for all samples, confirming the fractal self-similarity of the joint surfaces across the measured scale range. The slope of the PSD plot directly yielded the fractal dimension with values ranging from 1.26–1.50. Similarly, the BCM analysis produced fractal dimensions between 1.24 and 1.47, based on the scaling of the number of covering boxes. A direct comparison of the two methods showed excellent agreement, with a mean absolute deviation of 0.018 and linear regression$R^2=0.98$, thus validating the robustness and consistency of the fractal dimension as a quantitative roughness descriptor.

3.2. Fractional-Order Formulation of Damage Evolution

The fractional-order damage model, with derivative order α optimized to 0.82, accurately captured the progressive degradation of rock joints under shear loading. The stress–strain response simulated using this formulation matched the experimental data with a mean RMSE of 0.13 MPa. Figure 2 shows a numerical simulation of the fractional-order damage evolution in fractal rock joints under monotonic shear loading. Figure 2a shows the applied strain history, which simulated the standard direct shear test conditions. Figure 2b illustrates the damage variable $\omega \left(t\right)$ for $\alpha =0.75$ , demonstrating a smooth, non-linear increase that reflects progressive microcracking and weakening, which is a hallmark of post-peak softening in brittle materials. Figure 2c shows the sensitivity of the damage evolution to fractional order $\alpha$. Lower values (e.g., $\alpha =0.6$) produce overly diffused damage owing to enhanced memory effects, whereas higher values (e.g., $\alpha =0.9$) result in sharper localization, resembling classical local plasticity. The optimal value of $\alpha =0.7$ emerged as a material parameter directly linked to the degree of long-range correlation in microcrack development, validated by its ability to reproduce realistic post-peak behavior. The model successfully reproduced the initial nonlinearity, strain hardening, and post-peak softening behavior, demonstrating that the memory effect inherent in fractional derivatives is essential for modeling the history-dependent damage accumulation in geomaterials.

3.3. Existence and Uniqueness of the Solution (Banach Fixed-Point Theorem)

Applying the Banach fixed-point theorem to the integral form of the fractional constitutive equation yields a unique solution in the space of the continuous functions. Numerically, the contraction condition$‖{\rm T}\left(u_1\right)-{\rm T}\left(u_2\right)‖\le L‖u_1-u_2‖$ holds for all tested parameters with $L\le 0.85<1$. ensuring convergence of the iterative scheme and the good posedness of the model for numerical implementation. Figure 3 summarizes the results: Figure3a shows a realistic stress–strain curve, elastic loading, yield, and softening due to damage, in agreement with rock joint shear data; Figure3b shows fractional-order damage evolution ($\alpha =0.75$) capturing history-dependent microcrack growth; Figure3c shows linear plastic-strain accumulation typical of brittle materials under monotonic loading; Figure3d confirms rapid residual decay to machine precision, consistent with a contractive operator. These findings indicate a single, stable, and predictable response for admissible loading paths, and provide a theoretically grounded framework for geomechanical analysis beyond purely empirical models.

3.4. Lyapunov Stability of the Damage Evolution Law

Lyapunov functional $V\left(t\right)={‖D\left(t\right)‖}^2$was constructed for the fractional damage system. Using the properties of the Caputo derivative, it was shown that ${}_0{}^{c}D_t^{\alpha }V\left(t\right)\le 0$ for all $t>0$satisfies the criterion for asymptotic stability in fractional-order systems [41]. This result implies that the damage evolution process is stable and convergent; small disturbances in the initial damage state or loading path will not lead to unbounded or chaotic behavior, which is essential for reliable long-term geomechanical predictions. Figure4 shows the thermodynamic consistency and asymptotic stability of the proposed fractional-order damage model. Figure 4a shows monotonic damage growth without oscillations, whereas Figure 4b and Figure 4c confirm that the constructed Lyapunov functional $V\left(t\right)$strictly decreases over time with $V\left(t\right)\le 0$. This rigorous mathematical proof ensures that the model does not produce unphysical instabilities, even in the absence of the regularization terms typically required in local models. The results validate the use of fractional calculus as a physically grounded framework for modeling irreversible material degradation, thereby enabling reliable long-term simulations in complex engineering applications.

3.5. Enhancement of Shear Strength Parameters with Fractal Dimension

Quantitative influence of surface fractal dimension on key soil strength parameters. Specifically, Figure 5 a shows a nonlinear increase in cohesion with $D_S$, reaching a plateau, whereas Figure 5 b shows a linear increase in friction angle with $D_S$. As illustrated in Figure 5a, the effective cohesion ($c$) increases significantly with the fractal dimension ($D_S$), rising from 0.7 MPa at $D_S=1.24$to 3.5 MPa at $D_S=1.47$. This non-linear enhancement follows a power-law relationship, $c\left(D_S\right)=0.65+1.24{\left(D_S-1.20\right)}^{2.1}$, with a coefficient of determination $R^2=0.94$. Similarly, Figure 5b shows a clear linear trend for the internal friction angle ($\phi$), which increased from 31° to 49° over the same range as $D_S$, as described by $\phi \left(D_S\right)=30.5+1.20\left(D_S-1.20\right)$with $R^2=0.91$. These empirical correlations demonstrate that a higher surface complexity directly enhances both cohesive and frictional strength components.

3.6. Dynamic Evolution of Failure Surface Roughness

During the direct shear testing, the fractal dimension was monitored as a function of the shear displacement. Initially, the $D_f$ remained nearly constant during the elastic phase. As plastic deformation was initiated, $D_f$ increased by 6–9% owing to microcracking and asperity breakage, reaching a peak value. However, in the post-peak softening stage, intense abrasion and grinding caused the surface to smoothen, leading to a 4–7% reduction in $D_f$. This transient behavior, where $D_f$ first increases and then decreases, confirms its role as a dynamic state variable that evolves with the deformation history of the joint. Figure 6 demonstrates that the fracture surface roughness, quantified by fractal dimension $D_f$, evolves non-monotonically during shear, increases during propagation, peaks at peak strength, and decreases during localization. This dynamic behavior supports the treatment of $D_f$ as a state variable rather than a material constant. The 3D simulation confirmed the transition from roughening to smoothing, linking the microstructural evolution to macroscopic softening. These results provide a multiscale framework for modeling the strength degradation in geomaterials.

3.7. Fractal-Controlled Yield Surface in 3D Principal Stress Space

The modified Mohr–Coulomb criterion, incorporating fractal-dependent strength parameters, generates a yield surface in the 3D principal stress space that is substantially expanded compared with the classical model. For a joint with $D_S=1.40$, the fractal-enhanced surface predicts a uniaxial compressive strength 25% higher and a tensile strength 40% greater than the classical criterion using the same baseline $c$and $\phi$. This outward shift of the failure envelope reflects the additional energy required to overcome surface interlocking and asperity degradation, demonstrating the significant influence of geometric complexity on the overall rock strength. 

Figure 7 demonstrates how the microscale surface complexity governs the macroscopic yield behavior. The yield surface expanded outward with increasing $D_S$, reflecting the enhanced strength under multiaxial confinement. Crucially, the conical shape is preserved, ensuring compatibility with the standard plasticity theory. The volume of the elastic domain increased by approximately 20% as $D_S$ rises from 1.10 to 1.40, confirming that the surface roughness significantly amplified the load-bearing capacity. This visualization validates the proposed model as a physically grounded scale-adaptive extension of classical soil/rock mechanics.

4. Discussion

Introducing the fractal dimension into the Mohr–Coulomb failure criterion structure significantly advances quantitative knowledge regarding how surface roughness measures influence the mechanical properties of rock joints. Therefore, the assessment of the proposed modified criterion is presented in an unbiased manner in the following section, using SWOT analysis and comparison with traditional models. All subsequent claims in this study were justified and directly linked to the source of inquiry. The high point regarding the modified criterion is the high objectivity evaluation and the low levels of subjectivity. A critical shortcoming of the conventional JRC method is that it requires individual analysts to visually compare profiles with schematic drawings[35,36,37,38,39]. However, this method leads to user variability and low reproducibility, making it impossible to construct a base diagram. Conversely, the modified Mohr–Coulomb criterion offers a high-level overview of surface and strength parameter assessments. In this study, the box-counting analysis integrated into the new modified Mohr–Coulomb model displayed similarly high-objective measures. This level of data availability is an obvious advantage because it ensures or eliminates the engineer’s expert judgment profile when producing new parameters. This objectivity presents a significant asset, ensuring that the input to the new failure criterion is not skewed by an engineer’s judgment. Scale invariance and physical foundations. It is also essential to note the aforementioned issue of scale dependency of traditional JRC values, which do not have a true physical tie to the fundamental material properties. In contrast, the fractal dimension is truly scale-invariant, depicting the roughness for varying magnifications and sampling resolutions in a uniform manner. Moreover, our study discovered a clear cause-and-effect connection between this geometric characteristic and the intrinsic material characteristics.

The fractal dimension is an active state variable that governs the effective cohesion and internal friction angle through the mechanics shown in Figure 5. Thus, the fractal dimension ceases to be a mere descriptor and enters the realm of a fundamental geometric variable, introducing an immense level of physicality that is not present in empirical JRC. The modified model outperformed the traditional linear JRC adjustment by capturing the nonlinear saturation effect of strength enhancement at very high roughness levels, as shown in Figure 5a. Although the JRC-adjusted linear Mohr-Coulomb criterion often resulted in an unphysical strength decrease with increasing roughness, the modified criterion consistently incorporated a diminishing effect at high levels. Therefore, the modified model represents a more physical and versatile constitutive model that provides more accurate strength predictions over a broader range of rock joint conditions. Nevertheless, measurable geometric descriptions of the surface lead to linear results and significant weaknesses in terms of the high complexity and cost of the measurement, sensitivity to the measurement and processing, and an associated closed-form solution that might not be available in several cases. Although the findings showed a strong correlation for the materials tested, validation across a much more extensive range of conditions, such as cyclic loading, varying confining pressures, and different rock lithologies, is required. Despite their problems with objectivity and scale invariance, the classical JRC and Mohr–Coulomb have been used for almost a century and have been tested and validated on materials from all over the world.

This new approach may be theoretically and physically more accurate than the JRC method because it correlates the raw surface geometric data with mechanical constitutive parameters. However, a methodology is not necessarily superior when practicality and resource expenditures are considered. Hence, in many cases, the classical model is more than sufficient for engineering purposes because it is easier to implement and correlate with actual materials. The two major instances in which our model should be used instead are in highly detailed or research-oriented subsurface analyses, where surface roughness is of utmost importance, and in cases where litigation surrounding slope stability is important, a more detailed case could be made about the factors influencing the failure. Despite promising improvements in the predictive power for roughness-dominant failures, the proposed model has several current and perspective limitations. The current formulation assumes a deterministic correlation between Df and the strength parameters. However, because natural rock joints intrinsically exhibit spatial variability, it may be necessary to treat Df correctly as a stochastic quantity or upscale it to the level of field observations. In addition, it is necessary to confirm the prediction capabilities of the model for cyclic loading, fluid-saturated rock masses, and high-temperature environments. It may also be appropriate to extend the proposed model by introducing anisotropy into the fractal properties of the joint. Finally, the model should be validated for a larger number of rock types and joint conditions, and the formula should be further generalized for other rock types. Thus, even in the state of first approximation, the proposed model lays the groundwork for next-generation geomechanical models in which measurable geometric characteristics directly inform constitutive relationships to make more accurate and physically well-founded predictions of geotechnical stability than existing models.

5. Conclusions

This study presents a fractal-enhanced version of the Mohr–Coulomb failure criterion in which the surface fractal dimension directly appears as a state-dependent parameter governing the shear strength of rock discontinuities. By quantitatively establishing, with experimental certainty, the endogenous relations between the fractal dimension, effective cohesion, and internal friction angle, the model actively transmutes geometric roughness from a conceptually passive descriptor into a mechanically active variable. The application of dual box-counting and power spectral density fractal measuring methods provided both robust and objective input data, successfully eliminating the inherently subjective character of conventional assessments exemplified by the Joint Roughness Coefficient. Such complete and rigorous causal density determines the capacity of the framework to accurately capture history-dependent degradation and nonlocal behaviors through fractional calculus damage modeling, Lyapunov stability analyses, and existence-uniqueness proofs. The results demonstrate that the fractal-determined yield surface exhibits a substantial expansion with a higher fractal dimension, reducing the average mean squared error (AMSE)rate to less than 40% of classical formulations. The dynamic evolution of during shearing further confirmed its role as a state variable. This study establishes a scalable, physically grounded link between surface topography and constitutive behavior, offering a predictive foundation for advanced geomechanical analyses in rock engineering applications.