Overview of fracture definition – fractal geometry – fractal analyses

This work is an approach for studying tectonic events and their chronology during the Late Paleozoic period of Tassili-n-Ajjers in Algeria (Fig. 1), and underline the coexistence between the fractal dimension of the spatial faults distribution and the paleostrain tensors during the Hercynian tectonic event. In the published literature on this study area, only a few papers focused on this era, from the Early Carboniferous (the beginning of major Hercynian shortening) to the deposition of the first rocks of Triassic ages (358.9-251.9 Ma) (Latrèche, 1982; Boudjema, 1987; Haddoum et al., 2001; Perron et al., 2018; Chabane et al, 2019; Bosworth et al., 2020; Zazoun, 2020). The North America/Europe-Gondwana suture during the Upper Paleozoic was a principal geodynamic event in the Pangaea supercontinent formation and its history is recorded in a deformation zone stretching from Texas to Russia which includes the Ouachita, Alleghenian, Mauritanian, and Hercynian/Variscan (OAMH/V) fold belts (Riding, 1974; Hatcher, 2002; Franke et al., 2020; Nance et al., 2020). According Craig et al. (2008), across North Africa, the intensity of Hercynian deformation decreases from west to east. In Algeria, the Hercynian deformation is characterized by the reactivation of major sub-meridian Touareg basement fabrics inherited from the Pan-African Orogeny (900-600 Ma) (Caby et al., 1981). Thereby, on the uplifted highs; (e.g., Amguid–El Biod high, Hassi Messaoud ridge, Maouar high), the Cambrian–Ordovician sequence was affected by the Hercynian erosion and the anticlinal’s crests (e.g. Béchar, Ougarta, Timimoun, Sbaa, Bled El Mass, Reggane, Ahnet and Oued Mya areas) have been folded and eroded (Fig. 1A). (Boudjema, 1987; Craig et al. 2008; Galeazzi et al., 2010; Zieliʼnski,2011; English et al., 2016; Perron et al., 2018; Bosworth et al., 2020; Djezzar et al, 2020).

1.2. Overview of fracture definition?

Before studying the characterization of fractures, it is important to review some fundamental concepts. A natural fracture is a term adopted for designer for a planar macroscopic discontinuity, generated by brittle failure under an imposed stress field. Fractures in the rock are shown veins, joints and faults and they exist on a wide range of scales ranging from the microscopic scale (micrometer) to the macroscopic one (hundreds of kilometers). Based on our field observations and in accordance with the existing literature (Priest and Hudson, 1976, 1981; Fletcher and Pollard, 1981; Atkinson, 1987; Pollard and Aydin, 1988; Twiss and Moores, 1992; Aydin et al., 2006; Narr et al., 2006; Peacock et al., 2016), mineral-filled veins, joints and opened

fractures are referred to as mode-I (extension fractures). The faults correspond to mode-II fractures (Fig. 2A, B and C). The mode-III (tearing mode) refers to the mode where the relative movement is parallel to the fractures’ front. The Hybrid fractures combine both Mode-I and Mode-II or III.
Mode-IV (closing fractures) are anti-cracks and stylolites (Priest and Hudson, 1976, 1981; Fletcher and Pollard, 1981; Atkinson, 1987; Pollard and Aydin, 1988; Twiss and Moores, 1992; Aydin et al., 2006; Narr et al., 2006; Peacock et al., 2016).

The fractures description requires a concepts prioritization namely: from the most elementary (the individual fracture or fracture set), to intermediaries who are the fractures family then the fractures system, and up to the more global network. Fractures belong to a same system are the consequence application of a common stress or strain fields.

In outcrops, the joints appear as continuous and straight traces bearing no traces of displacement. According to the previous studies (Putot et al., 2001; Peacock and Mann 2005; Sorkharbi, 2014; Sanderson and Peacock, 2019), the joint networks can be morphologically divided into two types.
(1) In systematic joint networks: these are very regular fracture networks in terms of spacing and orientation (Fig. 2D).
(2) In joints with a very heterogeneous spatial distribution, as linear joint swarms (fracture corridors or fracture cluster), with strong sub-vertical development, on a more defined area, termed persistence (Fig. 2E).
In subsurface, fault networks and the natural fractures exert an important influence on fluids migration and the characteristics of natural fracture sets are related to tectonic events, diagenesis processes, and lithological variations (Laubach et al., 2010; Peacock and al., 2019).
In this study, we designate all the discontinuities under the term of fracture, which refers to a representation of discontinuities on two-dimensional map.

The fractures spatial distribution is often depicted as a fractal set, characterized by the fractal dimension (Bonnet et al., 2001). Using a filtered fault patterns, it is possible to establish the fractal dimension associated with the fracture patterns and the relationship between the field’s paleostress and the fractal dimension anisotropy (Pérez-López et al.,2005; Pérez-López and Paredes, 2006; Pérez-López et al., 2007). Indeed, in a fractured rock mass, the fractal dimension is a function of various parameters such as strength, tectonic history, and lithology. Fractal and multifractal methods have been widely used in geosciences and the latter provide precious information on the statistical and geometrical properties of geological parameters. The fractal geometry is applicable for the description of objects that show scaling behavior. For fractal geometry, the most important characteristic is the absence of a homogenization scale or representative elementary volume. Power law distributions are referred to as self-similar and their exponents are fractal dimensions. The notion of self-similarity is used for objects where parts of them look like the whole. However, power-law does not necessary imply self-similarity or fractality, however, some authors refer to power-law distributions as being fractal distribution and unlike mathematical fractals, geological phenomena exhibit fractal behavior within a limited scale range (Mandelbrot, 1982, 1985; Feder, 1988; Cowie et al., 1993; Bodin and Razack, 1999; Bour et al., 1997; Bour and Davy, 1999; Odling et al., 1999; Bonnet et al., 2001; Darcel et al., 2003; Bour et al., 2002; Rouai, 2006; Tran et al. 2006; Pérez-López et al., 2007; Xi et al., 2010; De Souza and Rostirolla, 2011; Kruh,l 2013). Scholz et al. (1993) found a linear scaling relationship between strain and length of faults, but they limited their model on faults that grow as isolated cracks.

1.3. Fracture properties affecting permeability

It is recognized that the permeability of the fracture network is governed by several factors, including orientation, aperture, height, length, and connectivity. In fact, the fractures are opened or sealed. However, the common cause of variability in permeability is the heterogeneous spatial arrangement (Nelson, 2001; Bisdom, 2016; Wang et al., 2019).

2. Objectives and Purpose

The purpose of this paper is to seek to establish the relationship between the fractal geometry of the spatial fault distribution and the strain field during the Hercynian tectonic deformation. In order to reach this goal: (1) we have mapped in the fracturing in area from lineament interpretation of aerial photographic dataset; (2) collected strike and dip fractures measurements from field work and applied fault population analyses to characterize the field’s paleostrain; (3) characterized fault traces in order to classify fractures into sets, based on tectonic events; (4) measured the fractal and the multifractal dimensions, and the fractal anisotropy of these fracture maps using box-counting method.
The questions highlighted in our case study are are: (1) are there any auto-similar properties for the analyzed fracture networks? (2) Is it possible to define a relation between tectonic event, strain field and fractal dimension?.

3. Geological Framework and Stratigraphy of the Study area

The area of study (Fig. 1) is located on the southeastern border of the Algerian Saharan platform commonly known as Illizi Basin, between the Berkine Basin to the north and the Touareg shield to the south. It is patterned in series of important geological entities, which are represented by en echelon folded structures (Figs. 1A and 1C). These structures are associated with wrenching structures such as (Assekaïfaf-El Adeb Larache-Tiguentourine-La Reculée, Tihemboka, Butte Noire, Dôme à Collénias, Edjeleh, Oued Oubarakat and Tan-E-Lak) (Fig. 1B). The zone of interest (Fig. 1C) covers the southern branch of the NE-SW trending Fadnoun fault. The latter is associated with a Panafrican basement fabric (the 08°30’E) (Fig. 1A). It is part of a large lineament network which extends to the north into the Illizi Basin and on where are built up many structures such as Assekaïfaf, El-Adeb-Larache, and Tiguentourine (Fig. 1B). A detailed stratigraphic study is not included in this paper as it focuses mainly on structural geology. However, a brief stratigraphic review is presented for a comprehension purpose (Fig. 3). The Paleozoic sediments lie unconformably over the Precambrian basement (Panafrican unconformity) and are unconformably overlain by the Mesozoic (Hercynian unconformity). Thus, the Paleozoic sedimentary strata consist almost entirely of siliciclastics with carbonates developed in the Carboniferous period. These are characteristic of a shield domain and stable tectonic conditions. The sedimentology of the Paleozoic series is characterized by braided fluvial systems building out braid-plain deltas into intracratonic Paleozoic Basins (Fekirine and Abdallah, 1998; Fabre, 2005)

Intrinsically, the 3-D strain data can be conveniently represented in a 2-D plot called Flinn diagram (Flinn, 1962) by using ratios of the principal strain axes (Log X/Y vs. Log X/Y; with X≥Y≥Z).
In order to quantify deformation in rocks, structural geologists use brittle structures markers (faults with slickenside lineation’s) to perform strain analysis. These last were analyzed at sixteen locations (Fig. 1C). These analyses helped to construct the paleostrain tensors axes (λ1≥λ2≥λ3, Finite strain ellipsoid axes; λ1: maximum direction of stretching, λ2: intermediate strain axis and λ3: maximum direction of shortening), supposing that the slip’s direction and orientation of on each fault reflect a single tectonic event. Fault populations were gathered on mechanical coherence criteria and field’s observations (eg. crosscutting relationships and the superposition of slickenside lineation on the fault plane). A numerical method (Sperner et al., 1993) was used to compute paleostrain configurations, defined by the orientations of the three strain tensors axes together, and to estimate the shape ratio (k) and the strain intensity (r) of the strain ellipsoid (Flinn, 1962). With k=1 (plane strain ellipsoid); k=0 (uniaxial oblate ellipsoid); k= (uniaxial prolate ellipsoid); 0<k< 1 (three axial oblate ellipsoid) and 1<k< (three axial prolate ellipsoid). Also, the strain values (in XZ, XY and YZ planes) have been used to evaluate the finite natural logarithmic strain ((Eq. (1)): Ramsay and Huber, 1983 and 1987

ε=(〖1/3)〗^(1/2) 〖[(〖ln⁡(R_XZ )〗^2+(〖ln⁡(R_YZ )〗^2+(〖ln⁡(R_Xy )〗^2]〗^(1/2) (1)

In total, 586 structural data sets (faults, bedding surface and fold axes) were collected during two field geology trip campaigns (2000 and 2002) in the Tassili-n-Ajjers area. 239 structural faults measurements from 16 field stations were obtained (Fig. 1A; Table 1). The paleostrain axes (λ1≥λ2≥λ3) reconstitution were established using the TectonicsFP © software (Ortner et al., 2002). The results of strain tensors parameters are reported in Tables 1 and 2.

4.1.2. Non-coaxial and coaxial shearing histories
A graphical method for evaluating kinematic patterns (Simple shear and pure shear deformation) at the crustal scale was developed by Gapais et al. (1987). This latest consists in plotting the shear zone data on a stereographic projection in a coordinate system defined by the finite strain axes (λ1≥ λ2≥ λ3). Following the definitions of Truesdell and Toupin (1960), the local geometry of a given shear zone is described by the shear direction (L), the normal (N) to the shearing plane (shear zone walls), and the normal (M) to the plane of shear (Gapais et al., 1987).

4.2. Structural Fault Map
The aerial photos analysis at scale of 1/80,000 resulted in the 2D lineament mapping of the Tassili-n-Ajjers area. This technique is used especially where intense fracturing, tabular structure and sub-vertical faults are underlined (Fig. 2D and E). The photogeology interpretation permits a suitable quantitative approach of fracture study.
The tabular geomorphology of the Tassili-n-Ajjers area and its lack of vegetated cover make it possible to avoid measurement problems.
The fracture trace maps files were prepared in two (2) pre-process stages:
Imported digital aerial photos into Adobe Illustrator™ software and drew fractures traces, which are underlined as ‘line’ or ‘polyline’ elements;
The layer is saved with the fracture traces as a Encapsulated PostScript file (.EPS).

4.3. Fractal Analyses
Since Mandelbrot’s work (1982, 1985) on the concepts of fractal geometry, several authors (Bonnet et al., 2001; Riley et al., 2011; Kruhl, 2013) have used the fractal analysis technique to characterize the two-dimensional geometry of fracture networks. In recent years, applications of fractals (Velandia and Bermúdez, 2018) and multifractals (Ord and Hobbs, 2019) analysis have been increasing in the earth sciences. However, the fractal dimension of the fracture network does not allow us to know the state of geometric connections between fractures (Bonneau, 2014).
In this paper, the fractal analyses approach consists to study (a) the fractal analysis (FA) of dynamic fault patterns; (b) the multifractal analysis (MFA) of field’s faults and (c): the analysis of the fractal anisotropy analysis (FAA).

4.3.1. Fractal Analysis
A fractal analysis (FA) was undertaken on the structural map to determine the faults’ spatial distribution and their relationship with tectonic events. Tectonic characterization of fracture traces can be carried out by classifying fractures into sets according to tectonic events. The filtering of different faults active under a particular strain field can be done by structural and tectonic analyses, such as faults’ sets analysis. By applying these techniques to filtering faults patterns initiated by each strain field, the fractal dimension obtained for these filtered maps (so-called dynamical maps) can provide the maximum roughness value for each fault set. Indeed, a Box counting method (BCM) (Roy et al., 2007) is a really suitable method for analyzing complex and isotropic 2D patterns, simple and non-centered, and which is not based on the central pattern symmetry (Kruhl, 2013). According to Paredes and Elorzaz (1999), the covering fractal dimension is determined by the classical box-counting technique through the study of the number N() of cells of size () that covers the set of fractures (F) when(Eq. (2)):

〖 N() 〗^(-D_b ) (2)

Thus, by reporting N() versus () in a logarithm plot, the fractal dimension D_b can be derived as the slope of a straight line. This method has been widely used to measure the fractal dimension of fracture networks (Bonnet et al., 2001).

4.3.2. Multifractal Analysis
The multifractal analysis concept was introduced in order to meet requirements that are not solved by the ordinary fractal (‘monofractal’) method. Thus, multifractal analyses (MFA) can provide more information on measurements of spatial objects than monofractal analyses. The traditional partitioning approach in the moment-based multifractal model is the box-counting method (BCM), in which the studied space is partitioned into non-overlapping boxes (Hentschel and Procaccia, 1983; Halsey et al., 1986; Evertsz and Mandelbrot, 1992; Cheng, 1999; Kruhl, 2013; Djezzar et al., 2020). For monofractals, scaling can be described by only one exponent (fractal dimension). This is not the case for natural fractal sets like fracture networks, which are multifractals. These objects can then be entirely described by a spectrum D(q) fractal exponents, the generalized dimensions, where the fractal dimension is D0 and the function D(q) is the multifractal spectrum (Turcotte, 1992). More references to the use of multifractal methods can be read on literature (Vignes-Adler et al.,1991; Davy et al.,1992; Sornette et al.,1993). According to Bonnet et al. (2013), the information obtained by the box-counting method, which only characterizes the scaling properties of the spatial occupancy of the fracture network, can be supplemented by the scaling properties of the fracture densities through the moments of order (q).
The system is first covered by a regular mesh of squares of side length () and the total length Li () of fractures in each square is measured. Then the probability pi () is defined as:

P_i ()=L_i ()/∑_1^n▒〖L_i () 〗 (3)

Where the sum is carried out over all boxes and gives simply the total cumulative length of all fractures. The moments of order q are then constructed (Eq. (4)):

M_q ()=∑_1^n▒[ 〖P_i ()]〗^q (4)

and should scale as M_q ()  ^( (q-1)Dq) , where the set {Dq for q= - to + } forms the “multifractal spectrum” of generalized fractal dimensions (Hentschel and Proccacia, 1983). Note that D_(q=0) = D by definition. The procedure of calculating the multifractal spectrum is conducted with q in the range from -20 to +20, with an increment of 1. The box dimension method and the multi-fractal analysis were achieved with the Software FracLac Version 2.0f ©2004 for Image J 1.40g (Wayne Rasband, National Institutes of Health, USA) developed by Karperien (2004).

4.3.3. Fractal Anisotropy Analysis
Several methods have been established to determine the anisotropy of structures in 2D. The structures deformation anisotropy is generally interpreted as strain. The analysis allows to construct a strain ellipse where the directional data are fitted to an elliptical distribution. (Volland and Kruhl, 2004; Kruhl, 2013). The technique proposed is simple. Firstly, they build a morpho-lineaments map from a digital elevation model (DEM) or aerial photos. Secondly, the fractal set is obtained by using Cantor’s Dust Method (Chilès, 1988; Velde et al. 1990) and compass-counting technique (Volland and Kruhl, 2004 Pérez-López et al.,2005; and Pérez-López and Paredes, 2006). Two parameters are obtained, the fractal dimension of spatial points distribution (D0) and the trend of the transect (δ). Plotting these parameters (D0, δ) in a polar graph allows determination of an ellipse, defined as the fractal ellipse of the lineaments’ spatial distribution.