Fractal Features of Fracture Networks and Key Attributes of Their Models

Abstract

This work is devoted to the modeling of fracture networks. The main attention is focused on the fractal features of the fracture systems in geological formations and reservoirs. Two new kinds of fracture network models are introduced. The first is based on the Bernoulli percolation of straight slots in regular lattices. The second explores the site percolation in scale-free networks embedded in the two- and three-dimensional lattices. The key attributes of the model fracture networks are sketched. Surprisingly, we found that the number of effective spatial degrees of freedom of the scale-free fracture network models is determined by the network embedding dimension and does not depend on the degree distribution. The effects of degree distribution on the other fractal features of the model fracture networks are scrutinized.

1. Introduction

Fracture systems are of paramount importance in hydrocarbon geology [1], hydrogeology [2], and geophysics [3]. The quantitative analysis of fractures and their patterns is essential for the evaluation of potential oil/gas production in reservoirs [4,5,6]. Fractures are habitually defined by their position, orientation, length, width, aperture, and roughness [7,8,9]. Additionally, the fracture system is often characterized by topological, morphological, kinematic, and hydraulic attributes [10,11,12]. Usually, these attributes are expected to be independent parameters, while the fractures form complex network systems [12,13,14,15,16,17]. In this regard, it was recognized that fracture networks often possess scale invariance [18,19,20,21] that allows for the upscaling of storage and transport properties of the studied fractured medium [22,23,24,25]. Accordingly, different models of fracture networks were developed using a fractal framework (see, for review, Refs. [26,27,28,29,30,31,32,33,34,35] and the references therein). In particular, the fractal geometry is well suited for describing multi-scale fracture networks under sparse data [35,36,37]. Indeed, the fractal approach allows the storage of the data relating to different scales of observation employing a relatively small amount of dimension numbers [38,39,40,41].

Although there is no canonical definition of fractals, the notion of a fractal network is commonly used in reference to a scale-invariant network whose fractal (e.g., self-similarity or box-counting) dimension D strictly exceeds its topological dimension d [41]. However, two dimension numbers (d and D) are often insufficient to properly characterize the striking properties of the fractal network. In fact, networks with the same similarity dimension $D>d$ can have very different topological, morphological, and metrological properties [41,42,43,44]. Hence, in order to unambiguously quantify the fractal features, one needs to employ more dimension numbers. In this context, we note that different fractal attributes have different effects on the network transport properties [45,46,47,48]. Therefore, a challenging problem is to define a suitable set of fractal attributes to account for the essential features of different fracture systems.

This work is devoted to the fractal geometry framework for the development of model fracture networks. The core aim is to establish a set of key attributes accounting for the essential fractal features of modeled fracture systems. The rest of the paper is organized as follows. In Section 2, the fractal features of fracture systems are discussed. The properties commonly measured in studies of fracture systems are briefly reviewed. The scaling features of fracture systems are highlighted. The fractal characteristics of fracture systems are outlined. The mapping from a three-dimensional fracture system into two- and three-dimensional fracture network is sketched in Section 3. A survey of the crucial features of fracture networks and corresponding fractal attributes is presented. A set of key fractal characteristics is suggested. In Section 4, we put forward two new kinds of fracture network models. The fractal features of these networks are scrutinized. The main findings are outlined in Section 5.

2. Fractal Features of Fracture Systems

The characterization of a fracture system in a permeable medium requires the analysis of the individual fractures (deformation bands, faults, cracks, stylolites, etc.) and their intersections. Indeed, individual cracks may or may not intersect. The intersecting fractures form fracture networks into three dimensions. Primary measures commonly used to characterize the fracture systems in soils and rocks are summarized in Table 1 (for more details, see Refs. [49,50,51,52,53,54,55] and the references therein).

A large number of experimental and theoretical studies suggest that there are similarity rules in the development of fracture systems [18,19,20,21,22,23,24,25]. This allows for the use of fractal geometry tools for the characterization and modeling of fracture systems. In particular, the scale-invariant fracture–pore space is commonly characterized by the fractal (box-counting) dimension D [38]. The fractal dimension also governs the sample size dependence of medium porosity, such that the overall porosity of the fractured porous medium scales as

$$\begin{array}{c}\hfill \varphi \propto {\left(\frac{{\mathsf{\ell }}_{max}}{{\mathsf{\ell }}_{min}}\right)}^{D-3},\end{array}$$

(1)

where $2<D<3$, while ${\mathsf{\ell }}_{min}$ and ${\mathsf{\ell }}_{max}$ are the characteristic sizes of the smallest and largest fractures (pores) in the system of size $L>{\mathsf{\ell }}_{max}$ [26].

Additionally, it was found that the fracture aperture can be linked to the fracture width by the empirical relation $A=c{W}^{\omega }$, where c is a fitting constant [19,50]. The deviations of scaling exponent $\omega$ from $\omega =1$ can be attributed to the self-affine roughness of fracture surfaces. The RMS roughness is defined as the standard deviation from the mean surface level ($\Delta h$). The self-affine roughness is usually characterized by the roughness (Hurst) exponent $\zeta$ defined via the scaling relation $\Delta h(\lambda \Delta x)={\lambda }^{\zeta }\Delta h(\Delta x)$, where $\Delta x$ is the window size on the mean plane and $\lambda >0$ [53]. The Hurst exponent is related to the box-counting dimension of the fracture surface as $D_B^{fs}=3-\zeta$ [11]. Furthermore, numerous experimental observations suggest that the length of fractures in the outcrops of fracture zones often obey a power-law distribution

$$\begin{array}{c}\hfill p\left(\mathsf{\ell }\right)\propto {\mathsf{\ell }}^{-\theta },\end{array}$$

(2)

where $\theta$ is the length distribution exponent [56]. Empirically, it was found that the values of the length distribution exponent are typically in the range of $1\le \theta \le 4$ [56,57,58,59].

The morphological features of fracture systems are frequently characterized by the lacunarity and the succolarity indexes [60,61,62]. Different fracture systems may have the same fractal dimension, but they can then be distinguished by lacunarity or succolarity.

Table 1. The common characteristics of fractures in soils, rocks, and geological formations.

Characteristic		Definition
Field observables	Fracture orientation	Spatial orientations in a sampling volume and on a sampling plane [49].
	Length (m)	Mean length of fracture traces on a sampling plane [12].
	Area (m${}^2$)	The area of the fracture plane [12].
	Volume (m${}^3$)	The volume of the fracture void [12].
	Aperture (m)	Distance between the two walls of a fracture [50].
	Spacing (m)	Spacing is defined as the averaged distance between the neighboring fractures in the fracture system [51].
	Intersections	Fractures’ intersections with a scanline and with a sampling area [52].
	Fracture surface roughness	Deviations of the fracture surface from the mean plane [53].
Intensive properties	Linear intensity (m${}^{-1}$)	Number of fractures per unit length [12].
	Areal intensity (m × m${}^{-2}$)	Fracture length per unit area [12].
	Volumetric intensity (m${}^2$ × m${}^{-3}$)	Fracture area per unit volume [12].
	Areal density (m${}^{-2}$)	Number of fractures per unit area [12].
	Volumetric density (m${}^{-3}$)	Number of fractures per unit volume [12].
	Porosity $0\le \varphi \le 1$	Porosity $\varphi$ is the ratio between the volumes of pore-fracture space to the volume of the sample [3].
Kinematic parameters	Displacement	The displacement of fracture walls against each other [12].
	Constrictivity factor $0<\delta \le 1$	The constrictivity factor is the arithmetic average of ratios between the areas of consecutive different cross-sections of the flow [48].
	Effective hydraulic aperture	Effective hydraulic fracture aperture is defined according to the cubic law [12].
	Filling	The fracture filling tell us whether a fracture acts as a conduit or prevents fluid flow [12].
	Formation factor $F={\tau }^2/\left(\delta \varphi \right)$	The formation factor can be determined as the ratio between the electrical resistivities of a fully saturated porous medium and the saturating electrolyte [54].
Topological features	Average degree of fracture system	Average degree of network $\langle k\rangle$ is equal to the ratio between numbers of fractures and intersections multiplied by 2 [55].
	Connectivity	Fracture system connectivity is commonly defined for a particular direction in terms of the relative fracture length projected into that direction [63]. If two fractures are directly connected to each other or there is a pathway from one to the other via other connected fractures, they have a connectivity indicator of one.

Table 1. The common characteristics of fractures in soils, rocks, and geological formations.

Characteristic		Definition
Field observables	Fracture orientation	Spatial orientations in a sampling volume and on a sampling plane [49].
	Length (m)	Mean length of fracture traces on a sampling plane [12].
	Area (m${}^2$)	The area of the fracture plane [12].
	Volume (m${}^3$)	The volume of the fracture void [12].
	Aperture (m)	Distance between the two walls of a fracture [50].
	Spacing (m)	Spacing is defined as the averaged distance between the neighboring fractures in the fracture system [51].
	Intersections	Fractures’ intersections with a scanline and with a sampling area [52].
	Fracture surface roughness	Deviations of the fracture surface from the mean plane [53].
Intensive properties	Linear intensity (m${}^{-1}$)	Number of fractures per unit length [12].
	Areal intensity (m × m${}^{-2}$)	Fracture length per unit area [12].
	Volumetric intensity (m${}^2$ × m${}^{-3}$)	Fracture area per unit volume [12].
	Areal density (m${}^{-2}$)	Number of fractures per unit area [12].
	Volumetric density (m${}^{-3}$)	Number of fractures per unit volume [12].
	Porosity $0\le \varphi \le 1$	Porosity $\varphi$ is the ratio between the volumes of pore-fracture space to the volume of the sample [3].
Kinematic parameters	Displacement	The displacement of fracture walls against each other [12].
	Constrictivity factor $0<\delta \le 1$	The constrictivity factor is the arithmetic average of ratios between the areas of consecutive different cross-sections of the flow [48].
	Effective hydraulic aperture	Effective hydraulic fracture aperture is defined according to the cubic law [12].
	Filling	The fracture filling tell us whether a fracture acts as a conduit or prevents fluid flow [12].
	Formation factor $F={\tau }^2/\left(\delta \varphi \right)$	The formation factor can be determined as the ratio between the electrical resistivities of a fully saturated porous medium and the saturating electrolyte [54].
Topological features	Average degree of fracture system	Average degree of network $\langle k\rangle$ is equal to the ratio between numbers of fractures and intersections multiplied by 2 [55].
	Connectivity	Fracture system connectivity is commonly defined for a particular direction in terms of the relative fracture length projected into that direction [63]. If two fractures are directly connected to each other or there is a pathway from one to the other via other connected fractures, they have a connectivity indicator of one.

Accordingly, the lacunarity index is frequently estimated using the gliding box algorithm for the calculation of the probability density functions of two- or three-dimensional images of the fracture system [64]. For a scale-invariant image, the probability density function $P(M,r)$ is a function of the box size r. The lacunarity index is commonly defined as

$$\begin{array}{c}\hfill \mathcal{L}\left(r\right)=\frac{Z_2\left(r\right)}{Z_1^2\left(r\right)},\end{array}$$

(3)

where $Z_q\left(r\right)={\sum }_M{M}^{q}P(M,r)$ are the moments of a density distribution function. For a two- or three-dimensional image of size L, the lacunarity index varies in the range between $\mathfrak{L}(r=L)=1$ and $\mathfrak{L}(r=1)=1/\varphi$, where $\varphi$ is the overall porosity. Accordingly, the normalized lacunarity index can be defined as

$$\begin{array}{c}\hfill \tilde{\mathcal{L}}\left(r\right)=\frac{\varphi \left[\mathcal{L}\right(r)-1]}{1-\varphi },\end{array}$$

(4)

such that $0\le \tilde{\mathcal{L}}\left(r\right)\le 1$ [61,65]. Alternative methods for characterizing the lacunarity of geophysical patterns were discussed in Ref. [66].

The succolarity characterizes a percolation capacity of the fracture system. Per definition, the succolarity determines the flow rate through the fracture network [61]. Accordingly, the succolarity for a given direction ($Su\left[dir\right]$) can be calculated as the normalized product between the fluid pressure and the area of the flooded space [64,65,67]. Explicit formulas for the estimation of succolarity from two- and three-dimensional images can be found in Refs. [61,64]. For anisotropic networks, the succolarities in different directions can be quite different. Accordingly, the two-dimensional anisotropy index is defined as the ratio of the succolarities in the horizontal and vertical directions

$$\begin{array}{c}\hfill \mathcal{A}=\frac{Su\left[BT\right]+Su\left[TB\right]}{Su\left[LR\right]+Su\left[RL\right]},\end{array}$$

(5)

where $dir=BT,TB,LR$, and $RL$ denote the bottom to top, top to bottom, left to right, and right to left directions, respectively, [61]. In a $3D$ succolarity calculation method, the virtual pressure fields are added to the $3D$ image in six directions [64].

A central issue for the transport properties of fractured porous media concerns the notion of the formation factor [68] introduced by Sundberg [69] and coined by Archie [70]. Conceptually, the formation factor accounts for the medium porosity ($\varphi$), the transport streamline constriction, and the tortuosity of transmission paths as follows $F={\tau }^2/\left(\varphi \delta \right)$, where the relative tortuosity of transmission paths $\tau$ accounts for the transmission paths lengthening, the network constrictivity factor $\delta$ accounts the variations in the streamline cross-section over the flow domain, while $\varphi$ restricts the area available for the mass transfer [71]. The path tortuosity is defined as the ratio between the actual path length to the shortest distance between the beginning and the end. The fractal geometry of the fracture system is reflected in the power-law dependence of the formation factor F on the fracture network porosity, widely known as Archie law [68]. Namely,

$$\begin{array}{c}\hfill F=c{\varphi }^{-m},\end{array}$$

(6)

where c is the fitting constant and m is the Archie exponent which varies in the range of $0.21\le m\le 3$ [72,73,74,75]. It has been argued that, for a scale-invariant fracture system, the Archie exponent can be expressed in terms of dimension numbers characterizing the fractal fracture network (as can be seen, for review, in Ref. [48] and references therein).

The topological features of a fracture system are associated with the system ramification, connectivity, and loopiness [41]. The ramification of fracture system can be quantified by the ramification exponent Q linked to the topological fractal dimension $D_{tB}$ [45]. The scaling behavior of the fracture system connectivity is characterized by the connectivity dimension $d_{\mathsf{\ell }}$ [76]. Furthermore, the connectivity dimension along with the fractal loopiness index ${\Lambda }_{loop}$ determines the numbers of the effective spatial ($n_{\gamma }$) and dynamical ($d_s$) degrees of freedom of random walker in the fracture system [77,78,79]. For readers’ convenience, the definitions of basic dimension numbers are summarized in the

Table 2. Dimension numbers and scaling exponents characterizing the fractal fracture networks.

Parameter (Symbol)			Definition
Dimension numbers	Box-counting dimension	$D_B$	Box-counting dimension is defined via the scaling relation $N_B\propto r^{D_B}$, where $N_B\left(r\right)$ is the number of n-dimensional boxes of size r needed to cover the fracture network, while $D_B\le n$ [76].
	Topological fractal dimension	$D_{tB}$	$D_{tB}=\min \{s:\exists X\subseteq \Omega \mathrm{such} \mathrm{that}{D}_B\left(X\right)\le s-1\mathrm{and}d(\Omega \setminus X)\le 0\}$ [41].
	Connectivity dimension	$d_{\mathsf{\ell }}$	$d_{\mathsf{\ell }}={lim}_{\mathsf{\ell }\to 0}ln{N}_{\mathsf{\ell }}/ln\mathsf{\ell }$, where $N_{\mathsf{\ell }}\left(\mathsf{\ell }\right)$ is the number of points connected with an arbitrary point inside of the $d_{\mathsf{\ell }}$-dimensional ball of radius ℓ around this point [76].
	Fractal dimension of the minimum path	$d_{min}$	The fractal dimension of the minimum path is defined via the scaling relation $\langle l_{min}\rangle \propto r^{d_{min}}$, where $l_{min}$ is the shortest distance between two randomly chosen points on the network, while r is the Euclidean distance between these points and $\langle \dots \rangle$ denotes the ensemble average [76].
	Fractal dimension of geodesic lines	$d_g$	The fractal dimension of geodesic lines is equal to $d_g=D_B/d_{\mathsf{\ell }}$ [79].
	Number of effective spatial degrees of freedom	$n_{\gamma }$	The number of effective spatial degrees of freedom is the number of independent directions in which a random walker can move without violating any constraint imposed on it by the network topology [77].
	Number of effective dynamical degrees of freedom	$d_s$	The number of effective dynamical degrees of freedom is equal to the spectral dimension, which is commonly defined as $d_s=-2{lim}_{\tau \to \infty }\{ln\langle P(r,\tau )\rangle /ln\tau \}$, where $P(\tau ,r)$ is the probability that a random walker on the network returns to its origin after $\tau$ steps, while $\langle \dots \rangle$ denotes the spatial average [79].
Scaling exponents	Crack length distribution exponent	$\theta$	The crack length distribution exponent is defined via Equation (2). Commonly, it varies in the range of $1\le \theta \le 4$ [56].
	Crack roughness (Hurst) exponent	$\zeta$	Roughness (Hurst) exponent is defined via the scaling behavior of the RMS roughness [53].
	Fracture aperture exponent	$\omega$	Fracture aperture exponent is defined by the scaling relation $\Lambda =c{W}^{\omega }$ and varies in the range of $0.4\le \omega \le 2.5$ [19,50].
	Network ramification exponent	$\mathcal{Q}$	The network ramification exponent is defined by Equation (7).
	Degree distribution exponent	$\chi$	The degree distribution exponent is defined via Equation (10).
	Fractal dimension of the random walk	$D_W$	The fractal dimension of the random walk is defined via the scaling behavior of the mean squared displacement of a random walker $\langle r^2\rangle \propto t^{2/2D_W}$ [77].
	Fractal dimension of the path tortuosity	$D_{\tau }$	The fractal dimension of the path tortuosity is defined via the scaling relation ${\tau }^2\left(\lambda L\right)\propto {\lambda }^{2(D_{\tau }-1)}{\tau }^2\left(L\right)$, where ${\tau }^2={\langle r^2\rangle }_L/L$ [48].
	Archie exponent	m	The Archie exponent is defined via Equation (6).

.

3. Survey of Fracture Network Modeling within Fractal Geometry Framework

Although, fracture systems are generally developed within a three-dimensional volume [80,81,82], the positions, orientations and relationships between individual fractures can either be mapped in three or in two dimensions for the convenience of the modeling procedure [83,84,85,86,87]. Accordingly, the fracture network can be treated as a system of branches and nodes embedded in the three- or the two-dimensional space. The nodes embody the tips intersections of fractures represented by branches. This allows the use of the graph theory for modeling fracture networks [15,87,88]. The percolation theory was also widely employed to model the fracture networks [89,90,91,92]. Another class of fractal networks used to model the fracture systems is constituted of the scale-free networks with a power-law degree distribution [93,94].

The topological features of the fracture network are the basis for understanding the correlations between the structure and transport properties of the fracture system. The key topological features are allied with the network ramification, loopiness, and connectivity [41]. The order of ramification of the network node j is equal to the minimum number of branches that are necessary to cut in order to separate an arbitrarily big set of nodes $C_j\subset \Omega$ connected to the node j. The order of network ramification is defined as $R=mi{n}_j\{R_j\left(C_j\right):C_j\subset \Omega \}$. Finitely ramified networks have a finite order R, whereas, in the case of infinitely ramified network, $R_j\propto L^{Q_i}$, where L is the size of the set $C_j\subset \Omega$, such that $R_j\to \infty$ as $L\to \infty$. Accordingly, the order of ramification can be characterized by the ramification exponent

$$Q=mi{n}_{j\in \Omega }\{Q_j\left(C_j\right):C_j\subset \Omega \},$$

(7)

such that, for the finitely ramified networks $Q=0$, whereas the infinitely ramified networks are characterized by $Q>0$ [45]. Furthermore, it was recognized that

$$Q=D_{tB}-1,$$

(8)

where $D_{tB}$ is the topological fractal dimension [41]. For example, the Sierpinski gaskets are finitely ramified ($D_{tB}=1$), whereas the standard Sierpinski carpets (frequently used to model fracture networks) are infinitely ramified and have $D_{tB}>1$ [45]. For the percolation clusters, the ramification exponent is equal to the fractal dimension of the red bonds $D_{rb}$ [95], such that the topological fractal dimension of the critical percolation cluster (CPC) is equal to

$$\begin{array}{c}\hfill D_{tB}=D_{rb}+1,\end{array}$$

(9)

whereas the topological fractal dimension of the backbone of the CPC is $D_{t{B}^{bb}}=1$, because the percolation backbone is finitely ramified per definition.

The fracture network connectivity is also characterized by the connectivity index $C_c$, the transitivity $C_t$, the betweenness centrality $C_{bc}$, and the clustering coefficient $C_{cl}$ (as can be seen in Table 3) together with the connectivity (chemical) dimension $d_{\mathsf{\ell }}$ (see Table 2) and the network degree distribution $p\left(k\right)$ [96,97,98]. The degree of a node is the number of other nodes (k) that it is linked to. Experimental observations reveal that the fracture networks frequently (but not always) obey a power-law cumulative distribution of the network degree

$$\begin{array}{c}\hfill P(>k)\propto k^{-\chi },\end{array}$$

(10)

where $\chi >1$ is the degree distribution exponent [93,94].

In this regard, we recall that the scale-free networks with $\chi >2$ can be embedded into a regular lattice via the minimization of the total length of the links in the system [98,99]. The formed clusters of successive chemical shells are found to be compact, while the dimension of the shortest path between any two sites is smaller than one. Specifically, it was established that a chemical distance in the scale-free network ℓ scales with the Euclidean distance in the embedding lattice l as $\mathsf{\ell }\propto l^{d_{min}}$, where

$$\begin{array}{c}\hfill d_{min}=\frac{d(\chi -2)}{d(\chi -1)-1}<1,\end{array}$$

(11)

is the fractal dimension of the minimum path in the scale-free network, while the network fractal dimension is equal to d [98]. Consequently, the connectivity dimension of the scale-free network embedded in a regular lattice exceeds the lattice dimension d.

Table 3. Attributes characterizing the fracture network connectivity.

Attribute	Definition
Connectivity index $0\le CI\le 1$	Connectivity index is the probability that two arbitrary points within the domain are connected [96].
Transitivity $0\le C\le 1$	The transitivity is the multiplied per 3 ratio between the numbers of triangles and connected triples in the network [97].
Betweenness centrality $0\le \mathcal{B}\le 1$	The betweenness centrality is defined as the ratio of the number of shortest paths passing through an edge to the total number of shortest paths between all possible pairs of vertices [8].
Clustering coefficient $0\le C_{cl}\le 1$	The clustering tells us how well a network is connected on a local neighbor-to-neighbor scale. The local clustering coefficient $C_{cl}^i$ is the ratio of the number of triangles involving vertex i to the number of connected triples having i as the central vertex. The global clustering coefficient $C_{cl}$ is the average over all the local clustering coefficients for each node [54].
Cyclic coefficient $0<\Theta <1$	The cyclic coefficient of a vertex is the average of the inverse of the sizes of the smallest cycles formed by vertex and its neighbors. The cyclic coefficient of a network is the average of the cyclic coefficient of all its vertices [42].
Fractal loopiness index $0\le \Lambda <1$	The fractal loopiness index is defined by Equation (12).

Table 3. Attributes characterizing the fracture network connectivity.

Attribute	Definition
Connectivity index $0\le CI\le 1$	Connectivity index is the probability that two arbitrary points within the domain are connected [96].
Transitivity $0\le C\le 1$	The transitivity is the multiplied per 3 ratio between the numbers of triangles and connected triples in the network [97].
Betweenness centrality $0\le \mathcal{B}\le 1$	The betweenness centrality is defined as the ratio of the number of shortest paths passing through an edge to the total number of shortest paths between all possible pairs of vertices [8].
Clustering coefficient $0\le C_{cl}\le 1$	The clustering tells us how well a network is connected on a local neighbor-to-neighbor scale. The local clustering coefficient $C_{cl}^i$ is the ratio of the number of triangles involving vertex i to the number of connected triples having i as the central vertex. The global clustering coefficient $C_{cl}$ is the average over all the local clustering coefficients for each node [54].
Cyclic coefficient $0<\Theta <1$	The cyclic coefficient of a vertex is the average of the inverse of the sizes of the smallest cycles formed by vertex and its neighbors. The cyclic coefficient of a network is the average of the cyclic coefficient of all its vertices [42].
Fractal loopiness index $0\le \Lambda <1$	The fractal loopiness index is defined by Equation (12).

The network connectivity and loopiness determine the numbers of effective spatial ($n_{\gamma }$) and dynamical ($d_s$) degrees of freedom on the fractal network [77]. The fractal loopiness of the fracture network can be characterized by the cyclic coefficient (see Table 3) or by fractal loopiness index defined as

$$\begin{array}{c}\hfill \Lambda =\frac{d_s}{n_{\gamma }}-\frac{1}{d_{\mathsf{\ell }}},\end{array}$$

(12)

such that the loopless networks have $\Lambda =0$, whereas networks with loops at all scales are characterized by $0<\Lambda \le 1-1/d_{\mathsf{\ell }}<1$ [41].

The percolation theory was used to study mass transport phenomena since its foundation. The main feature of the percolation processes is the existence of the minimum concentration of percolating elements for which a percolation cluster between the opposite sides of system. Accordingly, the percolation cluster can be used as a model of fracture network. The first attempts to model the fracture networks within the percolation theory framework were based on lattice percolation [100]. A simplest percolation process on regular lattices is a Bernoulli site percolation: each lattice site can either be open with probability p, or closed with the probability of $1-p$. The open sites represent pores. When the concentration of open sites exceeds the percolation threshold, there is a spanning cluster of pores (pore network) in which each pore is connected to at least one neighbor pore. The fractal properties of the critical percolation cluster (CPC) are known to be universal and determined by the lattice dimension [101,102,103]. Nonetheless, most of the percolation models of fracture networks are based on the continuum percolation with regard to the fracture distributions within the rocks (as can be seen in Refs. [104,105,106] and references therein). In these models, the medium permeability was calculated by triangulating each fracture and solving flow equations. In particular, the percolation theory framework was used to model 2D fracture planes with a power-law size distribution uniformly located in a 3D space [105].

The percolation processes on the scale-free networks were studied in Refs. [107,108,109]. It was found that the percolation transition exists if the degree distribution exponent is $\chi >2$. However, the scale-free networks with $2<\chi <3$ have a vanishing threshold. Conversely, the scale-free networks with $\chi >3$ undergo percolation transition at a finite threshold of dilution. Furthermore, for networks with a degree distribution exponent varying in the range of $3<\chi <4$, the critical exponents are found to be dependent on $\chi$ [107]. Specifically, it was found that the critical exponent $\beta$ characterizing the scaling behavior $P_{\infty }\propto {(p-p_c)}^{-\beta }$ of probability that an open site belongs to the spanning cluster is equal to

$$\begin{array}{c}\hfill \beta =\left\{\begin{array}{ccc}1/(\chi -3),& \mathrm{if}& 3<\chi <4\hfill \\ & & \\ 1,& \mathrm{if}& \chi >4\hfill \end{array}\right.\end{array}$$

(13)

independently of the embedding dimension d [102]. The correlation length critical exponent depends on the degree distribution exponent as

$$\begin{array}{c}\hfill {\nu }_{\mathsf{\ell }}=\left\{\begin{array}{ccc}(\chi -1)/(\chi -3),& \mathrm{if}& 3<\chi <4\hfill \\ & & \\ 3,& \mathrm{if}& \chi >4\hfill \end{array}\right.\end{array}$$

(14)

while the correlation length with respect to the intrinsic (chemical) metric scales as ${\xi }_1\propto {(p-p_c)}^{-{\nu }_{\mathsf{\ell }}}$ [102]. Other critical exponents and the fractal dimension of the critical percolation cluster can be established using the following scaling relations them

$$\begin{array}{c}\hfill D=d-\frac{\beta }{\nu }=\frac{1}{\nu \sigma }=\frac{d}{\tau -1}=d-\frac{d\beta }{2-\alpha },\end{array}$$

(15)

where the critical exponents $\alpha ,\sigma$, and $\tau$ control the scaling behaviors of the total number of finite clusters, the size of the largest cluster, and the finite cluster-size distribution, respectively, [101,102,103]. However, the fractal properties of percolation clusters in the scale-free networks were studied only above the critical embedding dimension $d_c\ge 4$ [108], whereas the fracture networks are embedded into three or two dimensions, and so the scale-free networks embedded in $d<d_c$ represent a special interest.

4. Modeling of Fracture Networks

As was pointed out in the previous section, the percolation theory framework has already been used to model the transport properties of the fracture systems. However, most of the fractal features of the fracture networks cannot be reproduced by the models based on the classical percolation on periodic lattices. In this work, we introduce two new kinds of models which allow accounting for some fractal features of real fracture systems. The first kind of model is based on the percolation of line slots which can have either a fixed size, or obey a power-law size distribution. The second kind of models explores the site percolation in scale-free networks embedded in the two- and three-dimensional lattices. The fractal features of introduced models are scrutinized.

4.1. Fracture Network Models Formed by Slot Percolation on Regular Lattices

In order to account for a finite length of individual cracks forming the fracture network, let us first consider the percolation of randomly distributed line slots in the regular lattice. A line slot is formed by n open adjacent sites aligned along one of the lattice axes. Thus, the line slot can be viewed as a straight crack of length $l=n>1$ and width $w=1$. Notice that, in contrast to the percolation of line segments studied in Ref. [110], the slots can intercross and overlap (see Figure 1), such that an open site can belong to two or more slots (cracks). The overlapping of slots aligned along the same axis form the crack of length $l>n$. Neighbor slots contacting along their direction form the crack of width $w\ge 2$. The contacts and crosses of slots aligned along different axes can be viewed as the vertices of the fracture network.

Numerical simulations in this work were performed on the square lattices of size $L^d$ with free boundary conditions. At each step of the simulation, the slot position and orientation were randomly chosen from uniform distributions, while the slot length is $n\ll L$. The line slot can cross the network boundary. In such a case, the effective slot length is reduced. When the concentration of open sites (voids) exceeds a critical value ($p_c$), the percolation cluster (fracture network) spans the lattice in one or more directions (see Figure 2a,b). In this regard, it should be noted that, due to the slot overlapping, intercrossing, and the boundary effect, the void concentration is always less than the ratio $Nn/L^d$, where N is the number of imposed slots.

The data of our numerical simulations reveal that, with the increase in slot length, the critical concentration of open sites decreases as $p_c\propto 1/n$. Consequently, the overall porosity of system with the fracture network (see, for example, Figure 2b) is considerably less than the overall porosity of system with the pore network modeled by the site percolation cluster (as shown in Figure 2a). We also noted that the results of numerical simulations are consistent with the assumption that the percolation of slots of fixed length belongs to the same class of universality as the ordinary lattice percolation. Therefore, the fracture network generated by the Bernoulli percolation of cracks (see Figure 2b) is characterized by the same dimension numbers as the site percolation cluster. Specifically, the fractal dimension obtained by the box-counting method is close to the expected value $D=1.89$. However, the normalized lacunarity index defined by Equation (4) and loopiness defined by Equation (12) of the fracture network differ from those for the pore network. Furthermore, we noted that the fracture network anisotropy defined by Equation (5) can be modified by the changing of the ratio between the numbers of slots aligned in different directions. The quantitative results for the crack percolation model will be reported after completing the comprehensive studies based on Monte Carlo simulations.

In order to account for the power-law distribution of the crack lengths in the fracture system, the length of randomly distributed line slots should be sampled from the power-law distribution

$$\begin{array}{c}\hfill p\left(n\right)=\frac{(1-{\theta }_0)n^{-{\theta }_0}}{n_{max}^{-{\theta }_0}-n_{min}^{-{\theta }_0}},\end{array}$$

(16)

with a suitable slot length distribution exponent ${\theta }_0$, while $n_{min}$ and $n_{max}$ are the minimum and maximum lengths of straight slots, respectively. The preliminary results of numerical simulations performed with this model reveal that the crack length distribution in the spanning cluster (see Figure 2c) also obeys the power-law distribution (2); however, the crack length distribution exponent characterizing the fracture network is $>{\theta }_0$, while the maximum crack length $l_{max}\to L\gg n_{max}$ as the system size $L\to \infty$. Furthermore, the preliminary results suggest that the fractal dimension of the minimum path in the fracture network obeying the power-law distribution depends on the slot length distribution exponent. Consequently, D, $D_{tH}$, and $D_W$ are dependent on ${\theta }_0$. The normalized lacunarity index and anisotropy index of the fracture network are functions of ratio $n_{max}/n_{min}$ and ${\theta }_0$. The detailed results of completed numerical studies will be published elsewhere.

4.2. Fracture Network Models Based on Percolation in Scale-Free Lattices

A fracture network with a power-law degree distribution can be modeled via site percolation in a scale-free network embedded in d-dimensional lattice. Accordingly, in this work, we explore the scale-free networks embedded in two and three-dimensional lattices, while the degree distribution exponent is varied in the range of $3<\chi <4$. The percolation transition in the scale-free network embedded in the d-dimensional lattice is characterized by critical exponents $\beta$ and ${\upsilon }_{\mathsf{\ell }}$ which are independent of d (Equations (11) and (14)). Conversely, from Equations (11), (14) and (15) it follows that the critical exponent characterizing the correlation length with respect to the Euclidean metric in the embedding Euclidean space depends on the lattice and depends on the lattice dimension d as

$$\begin{array}{c}\hfill \nu =\frac{{\nu }_{\mathsf{\ell }}}{d_{min}}=\frac{d{(\chi -1)}^2-\chi +1}{d(\chi -2)(\chi -3)},\end{array}$$

(17)

where $3<\chi <4$. The critical percolation cluster formed on the scale-free network embedded in the d-dimensional lattice represents a model scale-free fracture network. It is characterized by the same degree distribution exponent $\chi$ as the whole network. Furthermore, we noted that the fractal dimension of the red bonds in the critical percolation cluster is equal to $D_{rb}=1/\nu$. Accordingly, from Equation (9), it follows that the topological fractal dimension of the scale-free fracture network is equal to

$$\begin{array}{c}\hfill D_{tB}=1+\frac{1}{\nu },\end{array}$$

(18)

where $\nu$ is given by Equation (17). The graphs of $D_{tB}$ versus $\chi$ for $d=2$ and $d=3$ are shown in Figure 3a. It contrast to the topological fractal dimension, the topological connectivity dimension $d_{t\mathsf{\ell }}=1+1/{\nu }_l$ is independent of the lattice dimension d and is equal to

$$\begin{array}{c}\hfill d_{t\mathsf{\ell }}=\frac{2(\chi -2)}{\chi -1}.\end{array}$$

(19)

The graph of $d_{t\mathsf{\ell }}$ versus $\chi$ is shown in Figure 3b. The general property that the topological connectivity dimension can be either equal to or less than the topological dimension of the fractal network implies that the scale-free fracture network models have the topological dimension $d^{CPC}=1$ in two as well as in three dimensions.

The fractal dimension of the critical percolation cluster (representing the fracture network model) is determined by the first scaling relation in Equation (15). Accordingly, from Equations (13) and (17), the fractal dimension of the scale-free fracture network model is equal to

$$\begin{array}{c}\hfill D=\frac{d-d(\chi -2)}{d{(\chi -1)}^2-(\chi -1)},\end{array}$$

(20)

Figure 4a shows the graphs of connectivity dimension $d_{\mathsf{\ell }}=D/d_{min}$ for the scale-free fracture network model as the functions of $\chi$ for $d=2$ and $d=3$.

The transport properties of the fractal fracture network are strongly dependent on the numbers of effective spatial and dynamical degrees of freedom of walkers in the network [45,46,47,48]. The number of effective dynamical degrees of freedom is equal to the network spectral dimension. The admissible values of the spectral dimension are bound in the range of $2D/(D+d_{min})\le d_s\le D$ [41]. However, explicit expressions for the network spectral dimension are known only for a few kinds of fractal networks. However, it was argued [48] that the spectral dimension of the fractal network can be approximately estimated using the following empirical relation

$$\begin{array}{c}\hfill d_s=\sqrt{D_{tB}D},\end{array}$$

(21)

The graphs of $d_s$ versus $\chi$ calculated with the help of Equation (21) together with Equations (17), (18), and (20) are shown in Figure 4b. Notice that the spectral dimensions of the scale-free fracture networks are always less than 2, and so the fractal dimension of the random walk is $D_W>D$. This means that the random walks in fracture networks are recurrent. By following the arguments from Ref. [77], we find that the number of effective spatial degrees of freedom in the critical percolation cluster formed in the scale-free network embedded in d-dimensional lattice is equal to

$$\begin{array}{c}\hfill n_{\gamma }=2D-d_S,\end{array}$$

(22)

which differs from the relation $n_{\gamma }=2d_{\mathsf{\ell }}-d_s$ derived in [77] in view of the unusual condition $D<d_{\mathsf{\ell }}$. Surprisingly, we establish that the scale-free fracture network models are characterized by the universal numbers of effective spatial degrees of freedom determined by the network-embedding dimension. Specifically, we found that $n_{\gamma }=2.00\pm 0.04$ in $d=2$ and $n_{\gamma }=3.70\pm 0.05$ in $d=3$, independently of the degree distribution exponent varying in the range of $3<\chi \le 4$ (see Figure 4c). Conversely, the index of network loopiness defined by Equation (12) increases with the increase in degree distribution exponent $\chi$, as shown in Figure 4d.

5. Conclusions

In this work, we put forward two kinds of models for fractal fracture networks in two and three dimensions. The first kind of models is based on the Bernoulli percolation of line slots in regular lattices. The second one explores the site percolation in scale-free networks embedded in the two- and three-dimensional lattices. The inherent fractal features and key attributes of model fracture networks are outlined.

The numerical simulations of the slot percolation model reveal that the fracture network porosity decreases with increase in the slot length. Nonetheless, we found that the critical exponents and the fractal attributes of the critical percolation cluster are independent of the slot length. In order to model the power-law distribution of crack lengths in the fracture network, we perform the numerical simulations of the slot percolation model with the power-law distribution of the slot lengths. We found that the crack length distribution exponent characterizing the fracture network model is larger than the slot length distribution exponent. We also noted that the fractal dimension of the minimum path and other fractal attributes of the fracture network model are dependent on the slot length distribution exponent. A more detailed analysis of these observations will be reported elsewhere after completing a comprehensive set of Monte Carlo simulations.

On the other hand, we suggest that the scale-free fracture networks is modeled within a framework of site percolation on the scale-free networks embedded in two- and three-dimensional lattices. With this model, the effects of degree distribution on other fractal features of the model fracture networks are analytically revealed. In particular, we established the expressions for the fractal, the connectivity, the topological fractal, and the spectral dimensions in the terms of the degree distribution exponent and the embedding dimension. Thus, by varying the degree distribution exponent, one can model the scale-free fracture networks with a predesigned set of dimension numbers. Surprisingly, we found that the scale-free fracture network models are characterized by the universal numbers of effective spatial degrees of freedom which are determined by the dimension of the embedding lattice.

In summary, this work focused on the topological features of the fractal models of fracture networks. Our findings provide a novel insight into the modeling of fractal fracture networks. The ultimate aim of the fracture network modeling is to predict or reproduce the transport and storage properties of different fracture systems. Accordingly, the effects of the fractal features of the model networks on the transport phenomena and storage capacity will be the subject of forthcoming studies.