Drainage fracture networks in elastic solids with internal fluid generation

The first tiny bubble appeared in the gel after about 2 hours from the beginning of the experiment and 10–20 bubbles nucleated during the next 20 hours. Bubbles grew until they reached a diameter of approximately l₀ = 2.3 ± 0.2 mm and then changed into elongated cracks with sharp tips after 12 ± 3 min. However, for about 50% of the fractures, an initial bubble nucleation and growth stage was not observed with the image resolution used in these experiment (see movies in the Supplementary Material: crack_1_EPL_format.wmv, crack_4_EPL_format.wmv). Crack tips propagated with increasing velocity, at both ends of the fractures, unless they became blunted or closely approached the boundary of the gelatine layer or another crack. During the free propagation regime, before cracks started interacting with other cracks or with the boundary, the crack length evolved with time as

$$\begin{equation} l=l_0 \exp{(t/\tau)}, \end{equation} \tag{ 1 }$$

where t is the time since fracture nucleation. The increase in the timescale, τ with increasing gas production rate, γ, could be represented by τ = 0.54/γ − 0.20 ± 0.017 h. Each crack tip either propagated until it connected to another fracture or reached the gel boundary, or it slowed down and eventually halted inside the gel matrix, forming a dead end. When a fracture was not connected to other fractures or the gel boundary, its aperture a increased with time due to increasing gas pressure until a maximum aperture of approximately 3 mm was reached along most of the fracture length. When two cracks intersected their apertures changed rapidly indicating pressure equilibration. When a crack reached the gel boundary, its aperture collapsed as the gas drained from the system in less than one second. Once a crack formed, it continued to serve as an intermittent drainage pathway that closed and opened. New fractures nucleated and propagated until a steady-state fracture pattern stabilized, see fig. 1A.

To analyze the development of the fracture pattern, the images of fractured gelatine were processed such that all pixels that were once covered by a fracture were preserved in accumulated images that represent the development of the fracture network. The final accumulated fracture pattern for experiment Exp₁ is displayed in fig. 1B.

In most cases, both crack tips propagated linearly, with about the same velocity, until they reached a distance of λ_e = 10 ± 5 mm from another crack or a distance of 20 ± 5 mm from the gel edge. In such instances they would slow down, speed up, or change direction, depending on the circumstances. We will refer to λ_e as the elastic interaction length.

Using digital image correlation [17], the motion of tracer particles in the gel and of the gel edges was tracked. When cracks open, new volume is generated in the system, but no systematic expansion of the gel in the plane of the Hele-Shaw cell was found. However, by using an LVDT displacement sensor (Omega Engineering), we found that the glass plates confining the gel moved in the direction perpendicular to the plane of the Hele-Shaw cell plane, allowing the gel to increase in thickness. The total deflection of the plates corresponded to a total volume increase of (1.9 ± 0.1) × 10⁻⁶ m³. Since the gel is essentially incompressible, this corresponds to the volume of the gas trapped in the fractures and the largest total volume of the open fractures was estimated from image analysis to be (2.2 ± 0.1) × 10⁻⁶ m³.

To characterize the mechanisms of crack formation and growth we analyzed to what degree new fractures nucleate either in the middle of unfractured domains, near existing fractures, or randomly. When a crack formed it was assigned a number indicating the order in which it was formed as shown for experiments Exp₁ and Exp₂ in fig. 2A and B. Distance maps as shown in fig. 3 were computed where the colour indicates the distance r(x,y) from the coordinate (x,y) to the closest fracture or free boundary. For each nucleation event i, the distance r_i, position, (x_i,y_i) and gradient ∇r(x_i,y_i) were determined just before nucleation occurred. On a line through (x_i,y_i) in the direction ∇r(x_i,y_i), the local maximum of r is denoted r_max,i. The normalized distance r_i/r_max,i = 0 or 1 correspond to nucleation at an existing fracture or in the middle of an unfractured domain. The cumulative distribution N(r/r_max > r_i/r_max,i), which is plotted in fig. 4, characterizes the degree of randomness in the locations of the nucleation sites and is fitted by N = (r_i/r_max,i)^μ, μ = 1.94. μ ⩽ 1 (the exact value depends on the shape of the unfractured area) corresponds to completely random nucleation and $\mu \rightarrow \infty$ corresponds to deterministic nucleation in the middle of an unfractured domain. The inset shows the cumulative distribution of r_i and demonstrates that there is a minimum distance between nucleation sites and pre-existing fractures.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The direction of fracture propagation from the nucleation point was measured and it was found that the distribution of angles between the fracture directions and the distance gradients, ∇r(x_i,y_i), was uniform between 0 and 90 degrees.

The steady-state fracture network shown in fig. 1B was skeletonized and all intersections of fractures (blue circles in fig. 5A) and fracture tips (red squares in fig. 5A) were identified and labeled as nodes. A fracture may be divided into segments consisting of the part of the fracture between two adjacent nodes or between a node and an end of the fracture and these segments are identified by different colours in fig. 5A, (the outer boundary is considered to be one continuous segment).

Zoom In Zoom Out Reset image size

At each node we measured the angles between intersecting segments. There were two distinct peaks: one in the range 180 ± 5° and the other in the range 90 ± 10°. By joining segments that intersect at angles about 180°, we reconstructed the continuous fractures which formed by propagation from a single nucleation point. This fracture pattern has a temporal order of nucleation (see numbering in fig. 2). In addition a hierarchical order can be defined such that if fracture C ends at fractures A and B its generation number is n_C = max(n_A,n_B) + 1, where n_A and n_B are the generation numbers of fractures A and B [6]. The draining perimeter is considered to be generation 0. In the experiments, up to five fracture generations (n = 5) were measured. In fig. 5B, the fractures of Exp₁ have been colour coded to indicate their fracture generation numbers. A fracture is "younger" than the youngest fracture it forms a junction with. One possible alternative generation numbering scheme emphasizing the drainage aspect of the network would be n'_C = min(n'_A,n'_B) + 1, which corresponds to numbering of river segments in a drainage basin. Using this scheme the highest drainage generation is n' = 3.

The production of CO₂ provides the energy required for fracturing the gel and the fracture network evolution depends on the coupling between CO₂ transport and elastic interactions in the matrix. The observation that fractures initiate as bubbles far away from existing boundaries suggest that crack nucleation is not determined by elastic interactions. We therefore propose the following diffusion-nucleation model:

• 1)  
  Most of the generated CO₂ is dissolved into the gel, which consists of 95% water. The solubility and diffusion coefficient of CO₂ in the gel are approximately the same as in pure water, c₀ = 1.8 g/l and D = 1.85· 10⁻⁹ m² s⁻¹ [18]. The dissolved CO₂ has negligible mechanical effects on the gel matrix.
• 2)  
  Bubbles nucleate when CO₂ generation leads to local supersaturation in the gel. Nucleation is heterogeneous and only a small supersaturation is required.
• 3)  
  As bubbles grow they begin interacting elastically with the gel matrix and a transition from bubbles to fractures occurs.
• 4)  
  The subsequent fracture propagation is driven by the gas pressure on the fracture walls which induces elastic stresses in the material.

After a fracture has nucleated it propagates when the gas pressure exceeds some critical value p_c. When fractures propagate slowly the gas pressure is constant and equal to p_c. Using the ideal gas approximation $p_c = \frac {n_{CO_2}RT}{V}=\mbox{const}$, the change in fracture volume is

$$\begin{equation} \frac{{\rm d}V}{{\rm d}t}=\frac{RT}{p_c}\frac{{\rm d}n_{CO_2}}{{\rm d}t}. \end{equation}\noindent \tag{ 2 }$$

The fracture aperture is limited to a = 3 mm due to adhesion of the gel to the walls of the Hele-Shaw cell and the volume of gas in the fracture is V ≈lah, where l is the fracture length and h = 3 mm is the distance between the glass plates. During propagation of a single fracture we assume a constant CO₂ flux, ϕ, over the fracture surfaces

$$\begin{equation} \phi = \frac{{\rm d}n_{CO_2}}{{\rm d}t}\frac{1}{2lh} = {\rm const}. \end{equation}\noindent \tag{ 3 }$$

The expression

$$\begin{equation} l=l_0\exp\left(\frac{2\phi RT t}{p_c a}\right), \end{equation}\noindent \tag{ 4 }$$

for the fracture length can be derived by combining eqs. (2) and (3) and integrating. Here, l₀ is the integration constant, which characterizes the critical length above which a gas bubble evolves into a fracture. This corresponds to eq. (1), which was used to fit the experimental data and comparison of eqs. (1) and (4) indicate that 1/τ∝ϕ∝γ.

Various factors control the final configuration of the fracture network. Since the gas production rate is virtually constant for a long period after the last fracture has stopped propagating we assume that the system reaches a quasi-steady state in which the fracture network drains the gel matrix sufficiently to hinder further fracturing. Using a diffusion-nucleation model, the characteristic size of unfractured domains, for which loss of CO₂ by diffusion is sufficient to prevent nucleation of new fractures can be estimated.

Consider a gel segment bounded by two parallel fractures separated by a distance 2δ. Assuming a constant CO₂ production rate, γ, the steady-state equation for the concentration field in the gel is

$$\begin{equation} D\nabla^2 c +\gamma = 0. \end{equation}\noindent \tag{ 5 }$$

We assume that the flux ϕ of CO₂ across a fracture surface is proportional to the difference between the concentration c of dissolved CO₂ at the surface and the concentration required for chemical equilibrium with the gas pressure p = k_Hc in the fracture, i.e.,

$$\begin{equation} \phi=D \frac{{\rm d}c}{{\rm d}x} |_{x=\pm \delta} = k\left(\frac{p}{k_H} - c(\pm \delta)\right), \end{equation}\noindent \tag{ 6 }$$

where k_H is the Henry's coefficient and k is a rate constant of CO₂ evaporation from the gel. The solution of (5) is

$$\begin{equation} c(x) = \frac{p}{k_H} + \frac{\gamma \delta}{k} + \frac{\gamma}{2D}(\delta^2 - x^2). \end{equation}\noindent \tag{ 7 }$$

We assume that to nucleate a bubble a critical supersaturation c_c = p_c/k_H is needed. Inserting c(x = 0) = c_c into eq. (7) we may solve for the critical distance δ_c and find two scaling regimes:

• diffusion limited ( $k\rightarrow \infty$), $\delta _c=\sqrt {\frac {2D}{\gamma }\frac {p_c-p}{k_H}}\propto \gamma ^{-1/2}$,
• evaporation limited ( $D\rightarrow \infty$), $\delta _c=\frac {k}{\gamma }\frac {p_c-p}{k_H}\propto \gamma ^{-1}$.

It follows that in the diffusion-limited case γ^1/2δ_c should have a constant value and γδ_c should have a constant value for the reaction-limited (evaporation-limited) case. In Exp₁ the mean production rate was γ = 1.6 ml/h and 2δ_c = 22.2 mm (γδ = 17.76, γ^1/2δ = 14.0) and in Exp₄ the mean production rate was γ = 6.2 ml/h and 2δ_c = 10.7 mm (γδ = 33.2, γ^1/2,δ = 13.3). This is consistent with the diffusion-limited model, but not with the reaction-limited model. The shape of the cumulative distribution function (for r_i/r_max,i) shown in fig. 4, with μ = 1.94, is also consistent with the diffusion-limited scenario, but not with the reaction-limited model since in the reaction-limited limit μ ⩽ 1.

The short range of the elastic interactions is a consequence of the adhesion between the gel and the essentially rigid walls of the Hele-Shaw cell. A three-dimensional solid full of gas-filled fractures will also have a limited elastic interaction length due to the "screening" by open fractures. The elastic interaction between fractures in the model system results in the opening and closing of fracture apertures at junctions, which act as valves to the gas flow. This dynamic process, causing intermittent gas release, will be described in detail in a future communication. The adhesion between the gel and the glass walls of the Hele-Shaw cell also imposes an effective confinement on the system similar to the effects of compressive elastic stress on a three-dimensional system. Under such confining conditions, neighbouring fractures in analogous three-dimension systems will interact in a similar manner —when the aperture of one fracture increases due to increased fluid pressure it will increase the compressive stress on neighbouring fractures and this will reduce their apertures.

It should be noted that although δ_c ∼ r_min ∼ λ_e ∼ 10 mm in experiments Exp₁ and Exp₂, the similarity between λ_e and the other lengths is fortuitous. The elastic interaction length λ_e does not depend on the gas production rate as do δ_c and r_min and elasticity does not control δ_c and r_min.

The direction of gas flow follows the direction of the pressure gradient, but because confinement induces elastic interactions between fractures and nearby fractures and junctions, which may act as valves when the gas pressure is reduced, the direction of the pressure gradient is not static. Because of the "valve closing" effect, fractures that are connected at one end to a fracture network that reaches the edge of the gel may propagate at the other end until that end also contacts a pre-existing fracture. The end result is a network with loops. A hierarchical-fracture network [6, 7] is ideally fully connected.

A simple model that captures some of the essential features of our experimental system and incorporates river-like and hierarchical-fracture–like networks as limiting cases was developed (see fig. 6). A lattice model of size L × L was used, where L is the size of the gelatine layer. Fractures nucleate and propagate sequentially and form a "final" pattern when the system is completely drained, i.e., when the cumulative fracture length in the pattern exceeds or equals the cumulative fracture length l_t. l_t/L is 8.5 in Exp₁ and 8.7 in Exp₂. The nucleation sites are chosen according to a probability distribution P(x,y) that depends on the distance from the site to the nearest fracture or boundary. We chose a power law distribution such that the probability of nucleating a fracture at site (x,y) is

$$\begin{equation} P(x,y) = r(x,y)^{\mu-1}/\sum_{x,y} r(x,y)^{\mu-1}, \end{equation}\noindent \tag{ 8 }$$

where r(x,y) is the distance from the nearest fracture or the boundary. This choice is justified by the experimental distribution (fig. 4), where the power law exponent was determined to be μ = 1.94 ± 0.06. The black dots in fig. 4 show the resulting cumulative distribution of r_i/r_max,i of nucleation points for 20 realizations of the model. This demonstrates that although the model was constructed with a different metric from that used to measure μ an exponent of μ = 1.94 makes the model nucleation distribution lie within the range of uncertainty of the experimental nucleation distribution.

Zoom In Zoom Out Reset image size

Once a nucleation site is selected, the fracture orientation is randomly selected as horizontal or vertical. The fracture propagates by the same distance in both directions until it reaches and forms a junction with another fracture or the system boundary. At this instance the fracture may either cease to grow because it is drained by the joining fracture (with probability 1 − ω), or it may continue to grow until the other end meets an existing fracture (with probability ω). ω = 1 corresponds to all fractures connected at both ends which is typical of a hierarchical-fracture pattern [6]. ω = 0 corresponds to all fractures connected at one end only which is typical of river networks. The probability ω of a network can be deduced from the ratio between the number of dead ends, n_d, and the number of junctions, n_j: $\frac {n_d}{n_j}=\frac {1-\omega }{1+\omega }$ and ω are 0.41 and 0.32 for Exp₁ and Exp₂.

The generation and exsolution of CO₂ in a gel-filled Hele-Shaw cell generated a new type of drainage fracture network with topological properties intermediate between those of directed river networks and hierarchical-fracture networks. The valve-like opening and closing of fracture apertures at fracture junctions plays an important role in the formation of this type of fracture network. Diffusion of CO₂ in the gel inhibits the nucleation of new fractures in the vicinity of pre-existing fractures and a simple random nucleation model that takes this and the probability that fracture propagation stops before both ends of the growing fracture reach other fractures, reproduces key characteristics of the experimental fracture network.