Fractal Characterization of Multimodal, Multiscale Images of Shale Rock Fracture Networks

Abstract

An array of multimodal and multiscale images of fractured shale rock samples and analogs was collected with the aim of improving the numerical representation of fracture networks. 2D and 3D reconstructions of fractures and matrices span from 10⁻⁶ to 10⁰ m. The origin of the fracture networks ranged from natural to thermal maturation to hydraulically induced maturation. Images were segmented to improve fracture identification. Then, the fractal dimension and length distribution of the fracture networks were calculated for each image dataset. The resulting network connectivity and scaling associations are discussed at length on the basis of scale, sample and order of magnitude. Fracture network origin plays an important role in the resulting fracture systems and their scaling. Rock analogs are also evaluated using these descriptive tools and are found to be faithful depictions of maturation-induced fracture networks.

1. Introduction

Interest in the study of unconventional rocks (shales, mudstones, siltstones, marls, etc.) has significantly increased in the last decade. Historically, shale formations were regarded as either source or seal rocks [1] while being overlooked for their energy resource potential. This drastically changed with the increase in the use of hydraulic fracturing and horizontal wells. The subsequent surge of natural gas production in the United States brought unconventional formations to the attention of the energy community. Significant research has been conducted since to study unconventional porous media and its transport properties [2,3,4,5,6,7,8,9,10].

Even though the interest in characterizing and modeling shale rock properties continues to be relevant, environmental concerns regarding the control of greenhouse gas emissions and storage have prompted research on the use of depleted shale reservoirs for CO₂ storage. CO₂ capture followed by immobilization within geological structures is an option with the potential to decrease carbon dioxide emissions [11]. Shale formations are sites of interest for the subsurface storage of CO₂ because of their gas storage capacity within their organic matter and pore space [12,13,14,15,16]. Shale formations can also act as seals for potential subsurface hydrogen storage in conventional reservoirs [17,18].

In unconventional rocks, fracture networks are presumed to constitute especially significant pathways for fluid transport. Fracture networks are considered to play a relevant role in the flow dynamics of shale host rock [19,20]. Network connectivity in low-permeability geological formations such as shale is, therefore, a critical feature controlling fluid transport [21]. Common sources of fracture networks within shale are the product of natural geological evolution, organic matter maturation, and hydraulic stimulation, among others [22,23]. Additionally, fracture networks are expected to contribute to shale rock’s total CO₂ storage capacity, in particular around areas of improved injectivity [16].

Characterizing fractured unconventional rock remains challenging. The spatial heterogeneity and complexity of fracture networks at all scales compound the difficulty of creating physically accurate fracture networks and the models of transport within those networks. The multiscale nature of fracture networks and their complex geometry make efforts to characterize them accurately and efficiently a difficult task in order to help define basic transport properties [24].

The concept of fractal geometry was first introduced by Mandelbrot [25]. It has since been explored in many fields of science. Fractal geometry helps describe objects that are difficult to characterize using Euclidean geometry at all scales. The probability distribution N of a self-similar variable a conforms to an inverse power-law distribution [26]. The transition between scales depends upon D ( fractal dimension) and the slope k of a power-law curve as [27]

$$N\left(a\right)=k{a}^{-D}$$

(1)

Fractal shapes are by definition scale invariant and self-similar, and they are characterized by their fractal dimension. In geosciences, the physics of fracturing has been proposed to be self-similar and to follow the rules of fractal behavior. Many studies have applied fractal geometry to characterize and describe porous media and fracture patterns in natural rocks, particularly in major fault systems [28,29,30,31,32,33,34,35].

Numerous studies at various scales have shown that the distribution of some fracture properties, such as length, often follow a power-law and as such, do not exhibit a characteristic length scale. Fractal geometry has been acknowledged as a tool to characterize fractured systems [32]. Interdependence between fracture length distribution l and the fractal dimension often observed in most fracture networks may be central to determining network connectivity [21]. Several studies have used a fracture network model that delineates connectivity based on the assumption that the fracture length distribution n(l) follows a power-law with an exponent $\alpha$ such that:

$$n\left(l\right)=l^{-\alpha }$$

(2)

Under this model [32], if $\alpha$ > 3 for 2D networks, then connectivity is controlled by the smallest fractures and the system responds to the basis of percolation theory. Whereas, if $\alpha$ < 2, the large fractures dominate connectivity and transport and if 2 < $\alpha$ < 3, both small and large fractures control the network’s connectivity [21,36]. When this model is combined with fractal density distribution, the following results emerge.

• Length distribution prevails over fractal geometry, and large fractures dominate global connectivity for:

  $$\alpha <D+1$$

  (3)
• Connectivity is controlled by fractures much smaller than the system size when:

  $$\alpha >D+1$$

  (4)
• The system is self-similar when the transition between regimes and connectivity properties are scale invariant as

  $$\alpha \approx D+1$$

  (5)

The objective of this study was to characterize fracture networks in unconventional rock samples and their analogs by calculating fractal dimensions and fracture length distributions. The dataset consists of fractured shale rock and rock analogs. The multidimensional, multimodal, and multiscale imagery was sourced from published and unpublished data at sample scales ranging from 10⁻⁶ to 10⁰ m. The fracture network characterization was made on the basis of fractal dimension and the fitting of fracture length distributions as power-law functions. Image processing and segmentation were performed on each image of the dataset. The relation between the calculated fractal dimension and the power-law exponent of the fracture length distribution was used to evaluate the connectivity and scaling behavior of each fracture network across six orders of magnitude of sample size length, according to the connectivity model described above. The results are discussed on the basis of image scale, dimension, and fracture network origin. Finally, the feasibility of rock analogs and microfluidic models of shale rock fracture network systems are validated.

2. Materials and Methods

2.1. Sample Image Dataset

A collection of multidimensional, multimodal, and multiscale shale core sample and rock analog images with sample sizes ranging from microns to meters and resolution ranging from nanometers to microns was sourced from published [37,38,39] and our unpublished data. A visual summary of the data included in this study, sorted according to Euclidean (vertical axis) and length (horizontal axis) dimensions is shown in Figure 1. The images are otherwise not to scale. The image samples are grouped as grayscale microscopy, high-resolution photography, and tabular and graphic models according to imaging modality.

Table 1. Image data details.

Sample	Sample Size Length Units	Image Resolution	Euclidean Dimension	Image Modality	Fracture Origin	Data Provenance
B1	$\mathsf{\mu }$m	31.2 nm/vx	3	nano-CT	N	[40]
G1 *	cm	30.58 $\mathsf{\mu }$m/vx	3	micro-CT	HI	[37]
G2 *	mm	540 nm/px	2	SEM	HI	[37]
G3	mm	540 nm/px	2	SEM	HI	[37]
G4	mm	540 nm/px	2	SEM	HI	[37]
G5	mm	5.11 $\mathsf{\mu }$m/vx	3	micro-CT	HI	[37]
W1	m	0.3 × 0.3 × 2.4 mm/vx	3	tabular data	N, F	unpublished
W2	m	0.3 × 0.3 × 2.4 mm/vx	3	tabular data	N, F	unpublished
W3	m	0.3 × 0.3 × 2.4 mm/vx	3	tabular data	N, F	unpublished
W4	m	0.3 × 0.3 × 2.4 mm/vx	3	tabular data	N, F	unpublished
A1	cm	60.8 $\mathsf{\mu }$m/px	2	DSLR camera	N/A	[38]
A2	cm	60.8 $\mathsf{\mu }$m/px	2	DSLR camera	N/A	[38]
A3	cm	60.8 $\mathsf{\mu }$m/px	2	DSLR camera	N/A	[38]
A4	cm	60.8 $\mathsf{\mu }$m/px	2	DSLR camera	N/A	[38]
A5	cm	60.8 $\mathsf{\mu }$m/px	2	DSLR camera	N/A	[38]
A6	cm	60.8 $\mathsf{\mu }$m/px	2	DSLR camera	N/A	[38]
M1	mm	244.7 $\mathsf{\mu }$m/px	2	graphics	N/A	[39]
M2	mm	117.1 $\mathsf{\mu }$m/px	2	graphics	N/A	[39]

* images from the same sample.

presents a summary of the image data used in this study including scale, image resolution, imaging modality, and the origin of cracks and fractures: heat induced (HI)—via pyrolysis—or natural or fracked (N,F)). These distinctions do not technically apply to the rock analogs or models; therefore, those samples show an ‘N/A’ in that category.

Grayscale microscopy (G series): Scanning Electron Microscope (SEM) mosaics images and X-ray nano- and micro-tomography reconstructions of samples from the Green River formation were sourced from [37,40]. High-resolution photography (A series): These are images of rock analog materials at the meso-scale (lab-grade gelatin), acquired with a high-resolution DSLR camera in a setup and process described in [38]. Tabular and graphic (W series): Macro-scale fracture networks were generated from tabular fracture data reported on the Wolfcamp site in the Midland, TX Formation core samples described in [20], and illustrations of microfluidics designs (M series) at the micro-scale were sourced from [39]. While a number of fractures in the G series were heat induced, fractures in the W series originate from natural and hydraulic fracturing.

2.2. Sample Descriptions

This section offers further details on the relevant and available physical properties of the samples whose images were sourced for this study. This should help create a wider useful overview of the materials listed above. The information is broken down using the nomenclature in Table 1. A graphical summary of the collection of images is shown in Figure 1.

Sample B: The sample in the B category belongs to the Barnett shale play and was imaged on the X-ray nano-CT located at line 6-2c of the linear accelerator at SLAC National Accelerator Laboratory (SSRL). The sample is a 30 $\mathsf{\mu }$m diameter shale rock cylinder. Based on the reconstructed tomograms, it shows evidence of pyrite conglomerates as well as organic matter regions surrounded by an inorganic matrix and a fracture network. More details about this sample are found in [40].

G series. The samples grouped in the G series consist of four rock shale samples from the Uinta Basin, Green River Formation, considered thermally immature. Samples G1–G4 are described in great detail in [37]. Data from sample G5 is unpublished. Samples in the G series were selected for their varied hydrocarbon content and the result of heat-induced (HI) maturation. Details of their composition and other properties are in [37]. The 2D samples are physically rectangle-shaped while the 3D samples are cylindrical cores with diameters in the mm range. Images of these samples reveal laminations that alternate between large mineral and large organic matter contents. The organic matter content and its location within these samples is key to their reaction to HI maturation. In organic-rich layers, thermal maturation induces a volume increase driven by generated liquid and gas phases. Phase appearance in turn causes an increase in pressure that alters the amount and distribution of the strain energy of the rock. The dissipation of this energy is one of the causes of fracture formation. The migration of organic matter such as bitumen can also prompt the reactivation and expansion of local fracture networks [37]. The fracture networks developed as a result and the availability of the imaging made these samples relevant to our study.

W series: The samples grouped in the W series consist of 4-inch diameter Wolfcamp shale rock core samples extracted from the Hydraulic Fracture Test Site 1 (HFTS) in the Midland Basin, TX, USA. An extensive characterization of these samples is presented in [20]. Additionally, a non-published and quantitative description of the fractured cores—particularly in terms of depth, dip, and dip direction data—was utilized to calculate a reconstruction of the fracture networks as reported. This reconstruction was based on the geometry, dimensions, and orientation data reported for each fracture; therefore, the available data for series W were not instrument-acquired images but model reconstructions based on those parameters.

A series: The data in the A series consist of the images of shale rock analogs from [38]. In that study, an experimental setup consisting of a Hele-Shaw cell filled with a solution of gelatin, sugar, and yeast that developed a fracture network within the gelatin matrix caused by the consumption of sugar to produce carbon dioxide by the live yeast. This is an analog process for shale maturation. Lab-grade gelatins of varying material strength (S) were tested to assess the dependency of the fracture network characterization on the matrix’s geomechanical properties. As a result, a number of fracture network images were collected and six of them are used in this study. Three material strength levels were used and the samples shown here are grouped in pairs according to strength, such that S_A1 = S_A2; S_A3 = S_A4; and S_A5 = S_A6. Furthermore, they were grouped according to increasing material strength such that S_A5,A6 > S_A3,A4 > S_A1,A2.

M series: The data in the M series correspond to microfluidics devices created for the study presented in [39]. The objective of this study was to image solute transport within multiscale shale fracture networks. Two fracture networks were created for that study: (i) a microfracture network (here referred to as ‘M1’) with fracture apertures between 1 and 500 $\mathsf{\mu }$m, and a macrofracture network (here referred to as ‘M2’) with fracture apertures between 5 and 500 $\mathsf{\mu }$m. The geometry of the network and the angles between the fracture intersections is based on an actual fractured shale sample [39]. A graphic representation of these designs is shown in Figure 1.

2.3. Image-Based Fractal Dimension Estimation

2.3.1. Segmentation

The fracture networks of the porous media listed in Table 1 were pre-processed, binarized and then segmented prior to fractal analysis and fracture length distribution characterization. The methodology applied to each dataset group was specific to their modality. All image processing was performed in ThermoFisher’s Avizo^™.

Grayscale microscopy (B and G series): For the shale rock sample SEM images and X-ray tomography reconstructions, pre-processing included grayscale smoothing and de-noising filters in most cases. Segmentation was achieved through careful manual thresholding. Skeletonization (for 2D images) and fracture labeling were applied to the binary images in order to calculate the fracture length distributions.

High-resolution photography (A series): The images from these rock analog series were binarized after being pre-processed as described in [41]. Skeletonization and fracture labeling were applied in order to calculate fracture length distributions.

Tabular data and graphics (W and M series): The images for the W series were generated from tabular fracture descriptive data, i.e., depth, and dip, and dip direction angles. Images from microfluidics model designs from [39] were converted into 8-bit, binarized, skeletonized (2D images) and labeled in Avizo^™ in order to calculate fracture length distributions.

2.3.2. Box-Counting Method

Originally developed by the German mathematician Hermann Minkowski, the box-counting method estimates the fractal dimension of a space in an Euclidean context. In it, the binarized image of the fractal surface is completely covered by a number, N, of boxes of Euclidean size $\delta$_E. This process repeats itself over a range of box sizes $\delta$_E. If the system is fractal, the number N($\delta$_E) of such boxes is a power-law function of $\delta$_E. Fractal dimensions for 2D and 3D images were estimated using the ImageJ plugins FracLac and BoneJ, respectively. Fractal dimension values were double-checked with the ‘Fractal Dimension’ built-in function in Avizo^™.

2.4. Connectivity Properties

Fracture length distribution is an important parameter in determining the connectivity of fracture networks. When associated with fractal dimension, it may offer an important correlation that allows us to estimate connectivity regimes based on the scale of the dominant fractures of the respective system size [21]. In order to apply this criterion to our dataset, we calculated the fracture lengths for each of the segmented fracture networks by performing a labeling analysis. The sample Euclidean dimension dictates the particular process as follows.

Instrument-acquired images: Two-dimensional images are skeletonized prior to automatic labeling. Then, the ‘Label Analysis’ built-in function in Avizo^™ computes the fracture length for each connected component. Tomogram reconstructions can be automatically labeled in Avizo^™ and the ‘Label Analysis’ function computes the fracture length for each connected component.

Tabular data: Fracture length for the tabular data of the W series was obtained by calculating the hypotenuse of the dip angle, with the difference between the observed upper and lower fracture depth as the opposite leg of the triangle.

Length distribution calculations were performed with Mathworks’ Matlab^™ built-in functions. The power-law fitting was computed using supplementary Matlab^™ code provided by the authors of [42]. The power-law exponent $\alpha$ was calculated for each dataset.

If such a fit was possible, we evaluated inequalities (3)–(5) and used the connectivity characterization model proposed in [21] to designate a connectivity regime for each of our datasets and make inferences about the scaling of the fracture networks of the shale rocks systems presented in this study.

3. Results and Discussion

3.1. Fractal Dimension

Fractal dimension was calculated for each sample listed in Table 1. Fractal dimensions for three-dimensional datasets were converted into their two-dimensional equivalent using the following relation [21,32]:

$$D_{3D}=D_{2D}+1$$

(6)

in order to compare results across the full dataset as a group. Figure 2 shows data grouped by order of magnitude and summarizes the results as follows.

Image scale: The fracture network segmented for sample B1 has a fractal dimension equal to 1.73, and therefore it is not considered fractal because D_3D < 2. The average estimated fractal dimension ($\overline{D}$) for a subset of the G series samples (G2, G3, and G4) was 1.39, while $\overline{D}$_micro (samples G1 and G5) and $\overline{D}$_macro (samples W1–W4) have close averages values of 0.18 and 0.20, respectively. This result suggests a small reduction in fracture network complexity at scales larger than 10⁻⁴ and possibly a smaller frequency of fractures as well. These observations were qualitatively made during the segmentation stage of the samples as scales increased, as well as during the reconstruction stage of the macro-scale data.

Image dimension: In order to compare fractal dimension values between 2D and 3D, we used Equation (6) to downscale the three-dimensional values of D. Acknowledging that Equation (6) is an approximation [21], a decrease in fractal dimension is observed for 3D image data as compared to 2D. This decrease coincides with an increasing length scale, as detailed above.

Fracture network origin (rock shale samples only): On average, using pyrolysis to mature kerogen and organic matter (G series) seems to produce an only slightly larger fractal dimension on the studied samples than those representing the group of natural and hydraulically fractured rocks (W series). In our dataset, micro-scale samples (G series) were subject to significant heat leading to accelerated kerogen maturation. This process tends to originate fracture clusters near and on the high organic matter concentration zones, that are arranged in laminations for these samples. Therefore, the resulting fracture networks are at least partially correlated to the sample’s total organic content (TOC). These factors were extensively discussed in [37]. As the natural/hydraulically fractured samples are only represented at macro-scale, no definitive comparison can be made in this regard, other than the slight difference denoted above.

Rock and rock analogs. The rock analog samples included in the A series are described in detail in [38,41]. As previously mentioned, the main distinction between these samples is that the rock analog material varies in fracture toughness (S) such that S_A5,A6 > S_A3,A4 > S_A1,A2. Therefore, we subgrouped this set of samples into pairs of similar strength as described above. If we consider fractal dimension as a validation parameter for the rock analogs, the pair A5, A6 (strongest material) has the closest average fractal dimension value to the 10⁻⁵–10⁻⁴ m range (1.47 to 1.39, respectively). At the micro-scale (around 10⁻³ m), the A3–A4 pair (mid-strength) has the closest average to the rock shale samples of this study (1.30–1.18). The same result is obtained when the macro-scale samples (1.30 average for A3, A4 and 1.20 average for the W series) are compared. Fractal dimension is not the only validation parameter for rock analogs; however, it complements the connectivity assessment as shown in the following section. Fractal dimension helps determine whether an analog reasonably approximates some rock properties or transport behaviors.

There is not a very close match between the fractal dimensions of the microfluidics models (M series) from [39] and the rock samples analyzed in this study; however, following the above argument, their validity is warranted in other ways, as described in the experimental work of [39].

3.2. Connectivity

In this section, the discussion revolves around the interplay between the samples’ fracture network fractal geometry and the applicability of a power-law function to their fracture length distributions.

Table 2. Connectivity results.

Sample	Alpha ($\mathit{\alpha }$)	D	Connectivity Dominance
B1	3.45	1.73	small fractures
G1	2.01	2.11	transitional towards large fractures
G2	3.57	1.21	small fractures
G3	2.29	1.41	transitional toward large fractures
G4	2.64	1.54	self-similar, invariant
G5	2.74	2.25	transitional toward large fractures
W1	9.03	2.23	small fractures
W2	8.30	2.24	small fractures
W3	5.56	2.14	small fractures
W4	6.46	2.18	small fractures
A1	4.31	1.62	small fractures
A2	4.03	1.62	small fractures
A3	4.81	1.29	small fractures
A4	7.51	1.32	small fractures
A5	4.92	1.49	small fractures
A6	6.64	1.45	small fractures
M1	23.89	1.76	small fractures
M2	2.58	1.56	self-similar, invariant

summarizes the results of the main calculated parameters, $\alpha$ and D, as well as the inference made on each individual fracture network’s connectivity based on inequalities (3) through (5).

2D Shale rock images:The samples included in this group are G2, G3, and G4. The segmented fracture networks, histograms, and cumulative distribution function (cdf) curves for each sample are shown in Figure 3, along with dotted straight lines that represent the fit to power-law functions for each curve. As indicated in Table 1, the image modality is SEM for all samples. The segmented fracture network of sample G2 mostly consists of small evenly spaced fractures with an average size much smaller than the system size. As a result, sample G2 exhibits the smallest fractal dimension of the group (D = 1.21) and the largest value of $\alpha$ (3.57). Values of $\alpha$ > 3 suggest that connectivity is dominated by fractures of length much smaller than the system [21,24,32].

On the other hand, samples G3 (partial view) and G4 show a small number of large fractures that cut across at least half of the imaged area, along with more distinct fracture clusters with a wider range of fracture lengths, i.e., larger lacunarity, or the presence of significant fracture-free areas. The values of $\alpha$ are larger than 2 but less than 3. This result points to the equal importance of short and large fractures for network connectivity. The fracture system of sample G3 is an example of a fracture network where the loss of connectivity caused by fractal lacunarity is compensated by the presence and connectivity provided by one or more fractures of a length close to that of the system. This balance is only achieved if $\alpha$ < D + 1 [24]. In the case of sample G4, $\alpha$ ≈ D + 1, suggesting that short and long fractures have equal weight on sample connectivity.

In order to interpret these results, we must consider the samples’ compositional distribution, i.e., laminations and organic matter content distribution, as well as lamination orientation. In samples G2 and G4, laminations are horizontal, and in sample G3, laminations run vertical [37]. The lamination composition of these samples alternates between large and small organic matter content. It is also important to make the distinction that the segmented section of sample G2 belongs to a small organic matter content lamination. Because the fracture networks of these samples originated from heat-induced organic matter maturation, the laminations with larger organic matter content have larger fractures and fracture clustering. This explains the short length and evenly distributed fracturing of sample G2, and the clustered, large connecting fractures found for samples G3 and G4. The interpretation of their $\alpha$ and D values confirm these observations.

3D shale rock images—micro-scale: The samples included in this group are B1, G1, and G5. The segmented fracture networks, histograms, and cumulative distribution function (cdf) curves for each sample are shown in Figure 4, along with dotted straight lines that represent the fit to power-law functions for each curve. As indicated in Table 1, the image modality is nano-CT (also known as transmission X-ray microscope or TXM) for sample B1, and micro-CT for samples G1 and G5. The segmented fracture network of sample B1 mostly consists of small and jagged fractures of a narrow range of sizes, as its histogram shows. Similarly to sample G2, sample B1 exhibits a small fractal dimension (in fact, it is not a fractal system), and $\alpha$ > 3; therefore, fractures much smaller than the system size are expected to control connectivity. The fracture networks of samples G1 and G5 seem to have a dual fracture system. There are a small number of large fractures that span the diameter of the sample and match the samples’ large organic matter laminations, and a large network of small fractures, mostly located on the small organic matter content laminations. Both histograms confirm the observation and their $\alpha$ values indicate that the connectivity is dominated by both small and large fractures. Again, the fractal lacunarity of the samples is balanced out by the large fractures that span a dimension of the sample size (diameter) and the condition $\alpha$ < D + 1 is met.

3D shale rock images—macro-scale: The samples included in this group are W1, W2, W3, and W4. The segmented fracture networks, histograms, and cumulative distribution function (cdf) curves for each sample are shown in Figure 5, along with dotted straight lines that represent the fit to power-law functions for each curve. As mentioned above, the image data for these samples do not come from direct imaging but from data reconstruction based upon the recorded fracture geometry. As a result, the reconstructions used for this study carry some uncertainty associated with any missing or inaccurate tabular data. The average fracture for this dataset consists of incline planes that start and end at the edges of the core sample. The histograms for all samples show a relative homogeneity of fracture sizes. The power-law exponent $\alpha$ for all samples is large (≫3) and indicates that connectivity is dominated by fractures smaller than those of the system. In this case, given that $\alpha$ > D + 1, it is expected that the fracture network connectivity increases with decreasing scale [21].

Shale rock analog images: The samples included in this group are A1–A6. The segmented fracture networks, histograms, and cumulative distribution function (cdf) curves for each sample are shown in Figure 6, along with dotted straight lines that represent the fit to power-law functions for each curve. Gelatin of increasing material strength was used as a shale rock analog, as described in [38,41]. The resulting fracture networks exhibit behaviors correlating to the gelatin tensile strength and the histograms are skewed to the left, indicating a predominance of fractures smaller than the system size. The cdf curves and their power-law fit show values of $\alpha$ > 3, for all cases, suggesting alignment with results found for some rock shale samples. As with the macro shale rock series, $\alpha$ > D + 1, indicating that connectivity is sustained by fracture networks of length smaller than the system.

Microfluidic devices: The samples included in this group are M1 and M2. The segmented fracture networks, histograms, and cumulative distribution function (cdf) curves for each sample are shown in Figure 7, along with dotted straight lines that represent the fit to power-law functions for each curve. The design and use of these microfluidics devices was explained and detailed in [39]. The histograms show that the fracture networks of samples M1 and M2 are comprised of fractures of relative similar size. Sample M1 has a very large $\alpha$ value, mostly for larger fracture sizes, and along with $\alpha$ > D + 1 suggests that the smaller fractures dominate connectivity in this system. Sample M2 does not have an excellent power-law fit, but its $\alpha$ and calculated fractal dimension suggest that connectivity will likely be dominated by both short and long fractures (2 < $\alpha$ < 3) or be scale invariant ($\alpha$ ≈ D + 1).

Same sample: In Table 1, datasets G1 and G2 are marked with a * to denote that they belong to the same sample, imaged at the micro- and nano-scale, respectively. Figure 8 shows the cdf curves for samples G1 and G2, along with dotted straight lines that represent the fit to power-law functions for each curve. In this case, D_G2 is only slightly larger than D_{G1-2D equiv.} (1.21 and 1.11, respectively), but their $\alpha$ values are different (1.01 and 3.57, respectively). This difference may be partially due to connectivity distinctions between 2D and 3D networks, as well as probably due to the difference in organic content of each segmented sample area or volume, and its consequences, as previously described. Overall, it seems likely that the imaged sample area and image resolution are the main reasons the fracture length distributions diverge in this case.

Rock and rock analogs.Figure 9 shows a close up of cdf curves for the shale rock samples G2 and G5 shown in Figure 4 in addition to the curves representing the rock analog samples A1–A6 and M1 and M2 microfluidic devices, from Figure 6 and Figure 7, respectively. The curves were normalized by sample surface area. Solid and dotted trend lines were superimposed on some of these curves to characterize the connectivity and continuity between samples. In particular, the fracture network on sample G1 seems to correlate with the rock analogs A5 and A6 with the dotted trend line of $\alpha$ = 1.4; whilst sample G5 follows a connecting trend line to rock analogs A1 and A2, with the solid trend line of $\alpha$ = 1.8. These results validate gelatin as a rock analog for thermally induced fracture networks, consistent with the original premise [38,41].

The distribution curves for the microfluidic devices do not seem to match the shale rock curves by one or two orders of magnitude. An increase in the fractured area, while preserving the current fracture network lengths and designs, may be all the adjustment needed to shift these curves closer to the path of the shale rock curves, and thus enhance their modeling potential.

3.3. Multiscale Length Distribution

Shale rock only: Cumulative fracture length distribution curves are plotted in Figure 10 for all the shale rock fractal samples included in this study (the curve for sample B1 is therefore not shown). Curve values were normalized by sample surface area, and trend lines were traced in order to find the connectivity regime continuity throughout the scales [24,36] should they exist. The fractal fracture network scale range covered in this study spans the 10⁻⁶–10⁰ m interval. There is a lack of data for the length range between 10⁻² and 10⁻¹ m, that makes that scale range underrepresented. A data source for this scale could possibly be addressed with medical CT scanner images.

Two trend lines were traced over the curve groups. A solid line is shown over the curve for sample G2, and over the curves for samples W1–W4. The power-law exponent for this trend line is $\alpha$ = 3.6. A dotted trend line is shown over the curve group of samples G1, G3, G4, and G5, with a smaller slope equivalent to $\alpha$ = 1.6. We found that a single trend line would not be able to encompass the curves for all samples, and that samples could in fact be grouped by the fracture network origin, as denoted in Table 1. An exception to this grouping criterion is sample G1, that was exposed to heat-induced or artificial thermal maturation, but is not grouped with the other similar samples. Because this sample does not exhibit the same organic matter distribution as the other samples, it could be technically considered not part of the HI group. The cdf curve for sample G1 supports this claim, and the $\alpha$ = 3.6 trend line is a better fit than the $\alpha$ = 1.9 line.

These results suggest that the fracture networks within unconventional rocks are created in an important and distinctive factor with respect to their fracture network connectivity and scaling. For the natural or hydraulically fractured dataset group (samples G2 and W1 through W4), the power-law exponent $\alpha$ is always >3, indicating that connectivity is controlled by fractures much smaller than the system size. Furthermore, $\alpha$ > D + 1 suggests that the connectivity decreases with scale. On the other hand, the heat-induced organic matter maturation group (samples G1 and G3–G5), have $\alpha$ values such that 2 < $\alpha$ < 3, suggesting that both small and large fractures contribute to connectivity, although $\alpha$ < D + 1 points towards the dominance large fractures. The segmented images for these samples in Figure 3 and Figure 4 support these calculations, showing several fractures—of the same size as those of the system—that seem to run across the sample areas, connecting distant fracture clusters. These large fracture areas coincide with large organic matter content laminations.

4. Conclusions

Fractal dimension and the statistics of fracture length are geometrical information widely regarded as characteristic of most natural and induced fracture patterns. Knowledge of these key characteristics aids efforts to construct digital representations of fractured systems. In this study, we characterized a number of fracture networks imaged within shale rock samples and shale rock analogs by analyzing fractal dimensions and fitting fracture length distributions to power-law functions. We applied them as tools to characterize fracture networks originating from different processes: maturation due to the pyrolysis of organic matter components as well as natural and hydraulically induced fractures. The specific connectivity criteria applied here are explained at length in [21,24,32]. Specifically, we conclude the following:

• Throughout a sample size range spanning from 10⁻⁶ to 10⁰ m, most sample systems show a power-law exponent $\alpha$ > 3, suggesting that fractures one or more orders of magnitude smaller than the system size control connectivity and classic rules of percolation theory apply. These fractal fracture network samples also fall in the $\alpha$ > D + 1 category, indicating that connectivity decreases with system scale. These samples belong to the natural/hydraulically fractured group.
• Fractures in the samples matured using pyrolysis (G1, G3, G4, and G5) correlate with a different connectivity regime than the natural/hydraulic fractures. These fractures fit power-law exponents $\alpha$ such that 2 < $\alpha$ < 3, suggesting that a combination of small and large fractures interact to determine connectivity, but they also conform to $\alpha$ < D + 1, such that larger fractures prevail in system connectivity. This effect increases with system scale. Given the fracture length distribution shown on their histograms and their segmented images, it is probable that the larger fractures (of the same order of magnitude as system length) are more dominant than smaller fractures (one or more orders of magnitude smaller than the system).
• The distribution and abundance of organic matter of the shale samples has a direct and predominant effect on the distribution and clustering of the resulting fracture networks if the rock is (naturally or artificially) matured. In the samples included in this study, this effect resulted in large fractures located within the laminations of samples with a large amount of organic matter. These fractures dominate connectivity. In the case of hydraulic fracturing, the process is correlated to the mechanical features of the rock and the fracturing parameters that overall produces fracture networks with connectivity most likely dominated by smaller system size fractures.
• Gelatin is found to be an accurate rock analog for maturation-induced fractures as validated through multiscale connectivity trends.