Classification of fracture patterns by heterogeneity and topology

The geometry of a fracture pattern is controlled by the forces that act on the fracturing solid, and depends on rheomechanical properties and material heterogeneities. The forces that drive fracturing may be externally imposed and/or internally generated. For example, fracturing in rocks may be driven by externally imposed tectonic forces and body forces (gravitational forces), or by internal forces, often associated with volume alterations in the solid [1]. Volume alterations may stem from volume-changing chemical transformations, phase transitions or non-uniform thermal expansion or contraction due to heating or cooling.

In porous solids, deformation may be coupled with the fluid pressure field in the pore space, and poroelastic phenomena become important [2]. If the permeability is high, fluids are effectively transported by Darcian flow, driven by gradients in fluid pressure, which prevents localized pressure accumulation and large pressure gradients. If the permeability is low and the fluid pressure changes locally at a higher rate than fluid transport can compensate for, stresses may build up and cause fracturing. This in turn alters the permeability and may create pathways for fluids to migrate through or drain from the solid [3]. Cracking of mud and other unconsolidated porous materials during desiccation is a familiar example of fracturing driven by fluid withdrawal, and hydraulic fracturing is an example of fracturing caused by fluid injection. More than a million oil and gas wells have already been hydraulically fractured [4], and recent application of this technology to the recovery of shale oil and shale gas has become a topic of great economic importance and wide public and professional interest and concern [5]. Fracture patterns produced by fluid injection have been studied experimentally and numerically by several authors, for example [6–8]. Hydraulic fracturing may also be driven by internal fluid production associated with chemical transformations in the solid. Examples include release of fluids from hydrous minerals in subduction zones, which can trigger earthquakes [9], and methane released by fracturing during degassing of sediments, which may have implications for the Earth's climate [10]. Primary migration (the formation of hydrocarbon fluids in source rocks and the expulsion of the hydrocarbons into secondary migration pathways) is believed to be accelerated by pressure-induced fracturing when kerogen decomposes into lighter hydrocarbons [11–13].

The work described in this paper is concerned with fracturing in quasi–two-dimensional systems. Even in two-dimensional systems, fracture pattern geometry varies greatly, ranging from more-or-less polygonal or hierarchical patterns (osbserved in cracking mud, glazes on pottery or cracking paint films [1,14–16] and in two-dimensional simulations of rock weathering driven by solid volume-changing reactions [17]) to tree-like or branching patterns (observed in two-dimensional computer simulations and experiments on hydraulic fracturing [18]). A recent quasi–two-dimensional experimental model for fracturing during primary migration [19] generates fracture networks that are topologically intermediate between polygonal networks (with many loops) and trees (with many fractures terminating in dead ends). In these experiments, a layer of gelatin containing yeast and sugar is confined in a Hele-Shaw cell, and fracturing is driven by exsolution of dissolved CO₂ produced via fermentation. The intermediate topology is related to the complex dynamics of opening and closing of fractures and episodic expulsion of gas.

In the latter experiments CO₂ is generated more or less uniformly throughout the solid gelatin and diffuses down chemical potential gradients towards the boundaries of the Hele-Shaw cell and into fractures. Accordingly, fracturing initiates in the bulk, where CO₂ saturation is highest, and is inhibited near pre-existing fractures. As a result, there is no fracture branching, and fracture junctions only form when one fracture reaches another. In other systems, fractures nucleate at the external boundary or near existing fractures, for example when a brittle solid is cooled from the outside, when desiccation depends on transport of fluid to the exterior via open surfaces or when volume-changing chemical reactions localize near boundaries, where gradients in chemical potential are steepest [20–23]. In such cases, fractures penetrate from the outside of the solid towards the inside, and branching may occur. In nuclear fuel pellets both types of fracturing occur. The generation of fission gases generates fractures within the pellet [24], and the internal heating, which swells the inside of the pellet more than the outside, generates fractures that penetrate into the pellet from the outside [25]. In other systems, heterogeneity may play an important role, and can be a factor that leads to more randomly distributed nucleation.

In the next section we present a generic two-parameter model for simulating fracture patterns in two-dimensional systems. One parameter controls the ratio between the number of fracture junctions and the number of dead ends, quantifying the network topology. The other parameter controls the bias in nucleation position determined by the distance to the nearest fracture or external surface. This quantifies the heterogeneity of a fracture pattern, as nucleation close to existing surfaces leads to clustering of fractures. In the subsequent sections we demonstrate how model parameters may be recovered from the image of a generated pattern and show how the same scheme may be applied to quantify the geometry of real fracture networks.

The following generic model can be used to generate fracture patterns that vary in heterogeneity and topology. It is not constructed to mimic an actual fracturing process, although there may be instances where the driving mechanism can be related to the model parameters. Fractures are generated successively according to the following algorithm, illustrated in fig. 1:

• 1.  
  Select a nucleation site $(i,j)\in \{1,L\}\times\{1,L\}$ with probability

  $$\begin{equation} P_{i,j} = d_{i,j}^{\gamma}/Z, \end{equation} \tag{ 1 }$$

  where d_i,j is the distance to the nearest existing fracture or boundary and $Z = \sum_{i,j = 1}^{L} d_{i,j}^{\gamma}$;
• 2.  
  Select a random propagation direction;
• 3.  
  Propagate the new fracture in both directions

  • –  
    until each end meets an existing fracture or boundary (with probability ω) or;
  • –  
    until one of the ends meets an existing fracture or boundary (with probability $1-\omega$) (the fracture extends from the nucleation point by equal distances in both directions).

Zoom In Zoom Out Reset image size

The significance of the parameters ω and γ is illustrated in fig. 2, which shows examples of networks generated with this model. The parameter ω, which can be varied from 0 (trees) to 1 (polygons), is a measure of the network topology, and it can be related to the expected ratio of dead ends to junctions, R, by

$$\begin{equation} R = (1-\omega)/(1+\omega). \end{equation} \tag{ 2 }$$

The parameter ω quantifies the tendency of fractures to divide domains, which in turn may be related to the network's ability to relieve stresses or drain fluids. In the gelatin experiments of [19], for example, fractures tend to halt if they are able to expel gas through existing pathways. The parameter γ quantifies where new fractures nucleate relative to existing fractures, and may depend on loading, material heterogeneities and the mechanisms driving the fracturing process. For negative γ new fractures nucleate close to existing ones, so the process is surface driven and produces heterogeneity. For positive γ the process is bulk driven and leads to homogeneous patterns.

Zoom In Zoom Out Reset image size

Regularisation of (1) becomes a problem as the resolution L tends to infinity for negative γ. In this case Z diverges, indicating that fractures prefer to nucleate arbitrarily close to existing fractures. This is not a problem for our finite-resolution simulations, but the probability for nucleation close to a fracture is sensitive to resolution. In practice we believe that the case $\gamma\ll0$ is of limited interest, as nucleation near existing fractures dominates even for moderately low γ values. It should, however, be kept in mind that higher resolution is required to distinguish fracture patterns as γ becomes more negative.

Figure 2 illustrates that both γ and ω have visually distinct effects on the generated fracture patterns, suggesting that these parameters may be suitable quantities for characterizing real fracture networks. We introduced a special case of this model in [19] and were able to estimate model parameters from the time evolution of experimental networks, but in order to be more useful for characterizing real fracture patterns, we require that they can be estimated from a single image of a developed fracture network.

The quantity ω is readily found by counting the number of dead ends and junctions in the network and taking the ratio R. ω can then be recovered from eq. (2).

An approximate temporal hierarchy can be deduced from the intersections in a fracture network: If a fracture i terminates on a fracture j, it indicates that i existed before j. Fracture i is then called a direct descendant of j, and j is a direct ancestor of i. This observation is exploited in the algorithm for estimation of γ from images, which is described below and illustrated in fig. 3:

• 1.  
  Identify fractures and network connectivity:

  • (a)  
    Binarize and skeletonize the fracture network.
  • (b)  
    At each junction, determine which of the adjacent segments form part of the same fracture, and remove a pixel between continuous fractures and incident (descending) fractures. After this procedure, each fracture corresponds to a connected set of pixels (not necessarily on a straight line).
  • (c)  
    A fracture's ancestors are found by searching for fractures within two pixels from its endpoints.
• 2.  
  For each fracture $i\in \{1,2,\ldots,N\}$:

  • (a)  
    Make a binary image including all fractures but fracture i and its direct or indirect descendants;
  • (b)  
    Calculate the distance transform of this image, and rank points in the image according to their distance d_i,j from existing fractures;
  • (c)  
    Identify a set of putative nucleation points, $I_{\textit{nucl}}$:

    • If fracture i has one ancestor, let $I_{\textit{nucl}}$ be the center point;
    • Otherwise let $I_{\textit{nucl}}$ be the set of all pixels on fracture i;
  • (d)  
    For a list of trial values $\gamma^*$:

    • i.  
      Calculate the cumulative distribution $CP_i^{\gamma^*}\ (\text{rank})$ of nucleation probabilities by accumulating the probabilities $P_{i,j}= d_{i,j}^{\gamma^*}/Z$ for the ranked points.
    • ii.  
      Repeat i) but restrict the sum to putative nucleation sites to obtain $CP_i^{\gamma*}\ (\text{rank}|\text{I}_{\textit{nucl}})$.
    • iii.  
      Store $C_i^{\gamma^*}$, the value of $CP_i^{\gamma^*}\ (\text{rank}|\text{I}_{\textit{nucl}})$ at the point where $CP_i^{\gamma^*}\ (\text{rank})=0.5$.
• 3.  
  The most likely $\gamma^*$ value is the one that minimizes $\bigl|\sum_{i=1}^N (C_i^{\gamma*}-0.5)\bigr|$.

Zoom In Zoom Out Reset image size

Figure 4 shows the result of measuring γ from images of networks generated with the model. When the known temporal sequence is used, our algorithm successfully recovers the model parameters, with the exception of small γ values. When the approximate temporal hierarchy is used, there is a linear correspondence between the estimated and input values of γ, although there is a systematic overestimation of γ which increases with both ω and the input γ. The discrepancy is not surprising given that the distance maps used in the measurement algorithm assumes the existence of too many fractures at the time of nucleation. This shifts the distance map towards smaller distances and result in a bias towards higher $C_i^{\gamma^*}$ values in step 2(d)iii. of the measurement algorithm.

Zoom In Zoom Out Reset image size

Using least-squares fitting we determine that the discrepancy can be compensated for by transforming the measured γ value according to

$$\begin{equation} \gamma \rightarrow 0.7\gamma_e - 0.6\omega - 0.3. \end{equation} \tag{ 3 }$$

It may be possible to reduce the systematic errors by approximating the time evolution differently. Instead of keeping all fractures but a fracture's descendants and itself, we attempted to keep only direct ancestors, but this resulted in poorer estimates than the method proposed here. One could imagine more elaborate schemes, such as averaging over all allowed temporal sequences, but the number of possibilities grows rapidly with increasing system size, making it impractical. Tests using random temporal sequences resulted in apparently random γ values, showing again that the temporal sequence matters.

Figure 5 show three applications of the proposed classification scheme to natural fracture systems. The pattern in fig. 5(A) was produced during CO₂ expulsion from a layer of confined gelatin containing yeast that ferments sugar and generates CO₂, and it is taken from [19]. For this example we determine that $\omega \approx 0.6$ and $\gamma \approx 1.4$ using the procedure described above. The intermediate value of ω reflects the many dead ends in the network. They arise because fracturing is driven by gas pressure in the fractures, and once a pathway exists for the gas to escape through, further fracture propagation is inhibited. The positive γ is in accordance with the CO₂ transport in the gelatin being diffusion limited, as discussed in [19]; fractures nucleate far from existing fractures, where the CO₂ supersaturation is highest.

Zoom In Zoom Out Reset image size

Figure 5(B) shows ice fractures on Mars [26]. For this example we measure $\omega\approx0.5$ and $\gamma\approx0.3$. The pattern is slightly more heterogeneous than fig. 5(A), with more variation in domain sizes and fracture spacing. The patterns in fig. 5(A) and fig. 5(B) are topologically similar, but in some cases the fractures in fig. 5(B) are difficult to trace so this is subject to interpretation.

Figure 5(C) is taken from [27] and shows spheroidal weathering fractures observed along a road section at Nico Malan Pass, South Africa. For this example we measure $\omega\approx0.9$ and $\gamma\approx-0.1$. The pattern is hierarchical with sequentially subdivided domains, reflected in the high value obtained for ω. The sequential splitting of domains also accounts for the slightly negative γ value. For this example the fracture tracing is sometimes ambiguous, and some regions are finely fragmented. With better resolution one would therefore expect to measure a lower γ value.

In certain experiments the temporal evolution of a fracture system is recorded, and our scheme can be used to estimate γ based on actual distance maps at the time of nucleation. We did this for the gelatin experiments in in fig. 5(A) and obtained $\gamma = 2.3$, which is close to the value, $\gamma = 1.94$, reported by Kobchenko et al. [16]. The value we found by mapping the crack pattern onto the generic model and correcting with (5) was $\gamma = 1.4$, and even though both γ values suggest prounounced nucleation far from other cracks the numerical values differ. There may be several reasons for this difference: the statistical uncertainty stemming from the small number of fractures and the differences between the model and the real physical system. For example, some fractures in fig. 5(A) grow simultaneously, so the assumption of a strict sequential growth is not valid. In general we assume that the rate of nucleation is low compared to the propagation speed for the model to produced the correct dynamics. Even if this is not the case we can think of our scheme as a method to obtain quantitative measures that can be used to compare patterns. The details of the driving mechanism and evolution are factors that determine where the pattern is mapped in the model's parameter space.

We have presented a simple generic model for generating two-dimensional fracture networks, and the two model parameters can be used to classify fracture networks according to their spatial heterogeneity and topology. One of the parameters controls the bias in the position at which new fractures nucleate, away from or towards external surfaces or existing fractures, which adjusts the homogeneity of the fracture pattern. The other parameter controls the topology, defined by the ratio between the number of dead ends and the number of junctions in the network, and it can range from 0 (trees) to 1 (polygons). We devised a scheme for measuring the model parameters from images of fracture patterns and were able to recover the parameters from generated networks. We also classified three natural fracture systems to demonstrate the scheme.

Our work complements the vast amount of literature focusing on the scaling of unfractured domain sizes and fracture lengths, in the field, in laboratory experiments and in simple fracturing and fragmentation models [28–30]. The parameters ω and γ are, in principle, independent of the scale of the pattern. In practice, parameter estimates are better for larger systems with many fractures (due to better statistics), but more inaccurate if the unfractured domains are small compared to the lattice resolution.

The two parameters of our classification scheme are intuitively related to the appearance of a pattern, are relatively simple to measure and work for intermediate network topologies. However, we stress that the model is not intended to mimic the dynamics of a fracturing process, and it does not address all aspects of fracture network geometry. For example, the model does not respect the effects of elastic interactions, which would affect the direction of fracture propagation and intersection angles between fractures. To include this, one could bias the direction of fracture propagation depending on the direction of surrounding fractures. The assumption of random fracture directions is justified in systems where stress is localized by the attachment to a substrate, such as drying mud or paint films [14–16] or experiments in Hele-Shaw cells [18,19].

Certain features of fracture networks cannot be produced by the model, but can still be accommodated by the classification scheme: Firstly, the model prohibits fracture branching, but this can be handled in the classification scheme by grouping of branches. Secondly, the model evolves fractures sequentially and avoids situations where co-evolving fractures mutually descend from each other in the approximated temporal hierarchy. To resolve this, one may assume that fractures forming a loop in the temporal hierarchy nucleated at the same instant. An extension of the model to three dimensions is possible by representing fractures as planes propagating radially outward from nucleation sites. A decision would have to be made concerning, for example, what happens as fractures meet, but such details should not be a problem for the classification of natural networks, where ω may be defined as the ratio of the total length of the free fracture boundaries to the total length of intersection lines between pairs of fractures.

This work was supported by the Norwegian Research Council, through the Petromaks project and a Center of Excellence grant to the Physics of Geological Processes Center (PGP).