Fault damage zones have a higher upscaled permeability than the host rock because of a higher fracture intensity therein. Fracture distribution in the damage zone depends highly on the geometry of fault segments. However, precise images of architectural elements of large-scale faults at depth are difficult to obtain by seismic acquisition and imaging techniques. We present a numerical method that generates fault segments at multiple scales from an imprecise fault trace based on the fractal properties of these segments. The generated fault segments demonstrate hierarchical self-similar architecture, and their lengths follow approximately a lognormal distribution. These characteristics are similar to real fault segments observed in outcrops and seismic surveys. An algorithm that covers fault segments accurately with the minimum number of circles is proposed to calculate the fractal dimensions for both natural and computer-generated faults. The fractal dimensions of natural and generated fault segments are similar and range between 1.2 and 1.4.

Advancement in the understanding of slip along and fluid flow across faults requires a comprehensive knowledge of fault segments and their spatial distribution along the fault zone. Segmentation of different types of fault systems (strike-slip, normal, and reverse faults) plays an important role in finding the predominant locations and fractal dimension of fault splays [13]. These splays control properties of the damage zone around the parent fault and the geometrical and statistical properties of second-order fault networks around this fault [1, 46]. In general, the redistributions of slip and stress along fault systems and the resulting spatial distribution of damage zone are highly impacted by the fault architecture [1, 4, 7, 8]. Stepovers between segments control the development of fault core because of the high strain levels therein [1]. Fault core distribution is critical for fluid flow in the subsurface because, in general, damage zones have a higher upscaled permeability than the host rock, whereas fault cores usually have a lower permeability [1, 912].

Given that fault segments are components of the main fault, they are parallel to subparallel to the main fault and are described by echelon geometry. Faults are made up of many segments in both vertical and lateral dimensions [13]. The segments develop by partitioning a fault into numerous fault branches that appear to be unlinked and separated on a horizon of interest, whether on a seismic time slice or as fault traces on an outcrop [1421]. The discontinuous and segmented fault architecture has also been described by high-resolution seismic surveys [2228].

It is important to note the segmentation of a growing fault is opposed by the newer continuous links of the small segments that form a larger fault segment. This phenomenon is called the growth of faults by linkage and coalescence [7]. Interacting fault segments link as fault displacement grows [1416, 29]. Linkage occurs when the segments overlap during a substantial fault interaction through their associated stress fields [3032]. Such linkage and coalescence in a wide range of scales (mm-km) are observed from experiments [33] and outcrops [7, 34] and also simulated numerically [35, 36].

Through linkage and coalescence, small segments merge and form large segments. However, fault segments at a fine scale are still important, because damage zones show similar geometries across a wide range of scales (mm-km) [37, 38], and their spatial distribution firmly depends on the segment geometries. Only considering damage zones on a large scale underestimates complex structures on a fine scale and neglects their impact on the fluid flow [3941]. Precise images of the architectural elements of faults at depth are difficult to obtain by seismic acquisition and imaging techniques. Hence, seismic interpreters typically describe faults as single plane elements and fail to introduce the off-plane, closely spaced segments [1]. Therefore, the geometry of fault segments at depth is uncertain in general.

In this paper, we propose a method to generate detailed fault segments with only the trace length from an imprecise seismic map. Furthermore, we validate the correctness of generated fault segments from a fractal dimension perspective. We propose a modified algorithm to accurately calculate the fractal dimension of both natural fault segments observed from real seismic data or outcrops and generated fault segments. The generated fault segments are the premise for modelling associated damage zone around each fault segment and evaluating their impacts on the subsurface fluid flow.

This section introduces a method to generate fault segments from an imprecise seismic map and a modified algorithm to calculate their fractal dimensions.

2.1. A Method to Generate Fault Segments

Mandelbrot [42] proposed the concept of fractals, which can be used to characterize irregular sets, regardless of the scale at which these sets are examined. Fractal dimension is a measure of a set’s irregularity or complexity. Natural fault systems have fractal or multifractal characteristics [23, 4348]. Different types of faults are discontinuous and are composed of fault segments. Adjacent discrete segments step aside and overlap slightly to form an echelon fault geometry [22]. The echelon fault geometry can be found at different scales, shown in Figure 1. This phenomenon depicts the fractal characteristic or the hierarchical self-similar architecture of fault systems [1, 33]. Furthermore, from experiments [33] and outcrop observations [33, 37], such self-similarity shows a wide range of scales from millimeters to kilometers. There are always two to five segments in any given hierarchical rank, and three or four segments are most observed.

Figure 1

(a) Echelon breaks in the Imperial fault from [49]; (b) stepping arrangement of the San Andreas fault at the Gualala River from [50]; (c) the Calaveras and Hayward-Rodgers Creek fault zone simplified from [51].

Figure 1

(a) Echelon breaks in the Imperial fault from [49]; (b) stepping arrangement of the San Andreas fault at the Gualala River from [50]; (c) the Calaveras and Hayward-Rodgers Creek fault zone simplified from [51].

Using outcrop images, de Joussineau and Aydin [1] investigated segmentation along strike-slip faults and found correlations between different parameters of the fault segments. These parameters are illustrated in Figure 2. According to the study of de Joussineau and Aydin [1], there exist power law relationships between the
  • (1)

    Mean segment length (L¯se) and the maximum fault offset (S)

  • (2)

    Mean step length (L¯st) and the maximum fault offset (S)

  • (3)

    Step length x and step width y

Figure 2

Geometrical parameters of fault segments used in this study. “Subfault” denotes a younger fault generated from an older one.

Figure 2

Geometrical parameters of fault segments used in this study. “Subfault” denotes a younger fault generated from an older one.

Generation of fine-scale fault segments from an imprecise seismic map is an inverse problem. The fault segments generated from a highly underdetermined inverse problem are nonunique. Our method incorporates minimum information (the trace length from an imprecise seismic map) to generate a detailed structure of fault segments. If more information is available, such as the slip distribution, it can help to constrain the results and make them closer to reality. To implement the algorithm, we need to make the following assumptions:

  • (1)

    Fault segments are self-similar

  • (2)

    Subfault segments of a younger generation are independent of each other

  • (3)

    The relationships found in the study of de Joussineau and Aydin [1] are valid at each generation step of the subfault segments

If these assumptions are valid, the fault segments of a younger generation should be constrained by Equation (4).

where the length of the imprecise fault trace from a seismic map, L; the number of segments in each generation, nse; the mean segment length, L¯se=1/nsei=1nseli (li is the length of the subfault i); the mean step length, x¯=1/nse1i=1nse1xi (xi is the step length of the ith step); and the step width of the ith step, yi, are all defined in Figure 2. S is the maximum fault offset. We assume that the step length is proportional to its width with a factor of 2.7.

The fault trace length, L, in Equation (4) can be obtained from the fault map, and we can solve for three unknowns (S, x¯, and L¯se). We assume that there are three or four segments in each generation with an equal probability of 0.5, respectively. To make the generated fault segments closer to reality, we include both overlapping and underlapping between segments. The overlapping structure is dominant with a probability of 0.9. However, to obtain the length of each subfault li, we need to assume a relationship between the length of each subfault and the mean segment length. We assume that the subfault lengths are either equal, li=lj=L¯se, or that their length ratio is li/L¯se=rand0.8,1.2 (a random number uniformly distributed between 0.8 and 1.2). These assumptions are implemented for individual step length (xi) and the mean step length (x¯). The fault segments generated with both assumptions are shown in Figure 3.

Figure 3

Three generations of fault segments obtained from an imprecise fault trace. The black line segment denotes an imprecise fault trace taken from a seismic map; green, blue, and red line segments are the fault segments. In each generation, the subfaults are either (a) of equal lengths, li=lj=L¯se, or (b) of nonequal lengths, li/L¯se=rand0.8,1.2.

Figure 3

Three generations of fault segments obtained from an imprecise fault trace. The black line segment denotes an imprecise fault trace taken from a seismic map; green, blue, and red line segments are the fault segments. In each generation, the subfaults are either (a) of equal lengths, li=lj=L¯se, or (b) of nonequal lengths, li/L¯se=rand0.8,1.2.

To the best of our knowledge, there are no field data reporting equal lengths of fault segments. de Joussineau and Aydin [1], Davy [52] found that the fault segment lengths obey roughly a lognormal distribution. We assume that (i) for a given subfault generation, the length of a fault segment is independent of its neighbours, (ii) the set of subfault lengths for a younger generation is obtained from an older one through a recursive division, and (iii) the length of the older segment is approximately equal to the summation of each segment in the younger generation. The length of the nth generation subfault is then calculated in Equation (5).

where X1,X2,Xn are the independent random variables from a uniform probability distribution. If we take the logarithm on both sides, logLn is the summation of n independent random variables logXi. According to the Central Limit Theorem, logLn should follow a normal distribution and, consequently, the subfault lengths, Ln, should follow a lognormal distribution (example in Figure 4(a)). However, n should be large to make the Central Limit Theorem hold. Also, the standard deviation of Ln decreases with an increase of n. The lognormal distribution of segment lengths is found in one generation of fault segments. If the fracture lengths are measured in a target domain with many generations of fractures, it is more likely to follow a power-law distribution [53] because of the universal principle for fractality of solids [54].

Figure 4

Lengths of generated fault segments from a 10 km fault trace with the lognormal fit: (a) 5 generations of the fault segments without merging (2238 segments in total and L¯se=192 m); (b) 5 generations of the fault segments with merging (1709 segments in total and L¯se=215 m).

Figure 4

Lengths of generated fault segments from a 10 km fault trace with the lognormal fit: (a) 5 generations of the fault segments without merging (2238 segments in total and L¯se=192 m); (b) 5 generations of the fault segments with merging (1709 segments in total and L¯se=215 m).

The self-similar fault segments can have many subdivisions until reaching the scale of microfractures. Damage zones also show similar geometries across a wide range of scales (mm-km) [37, 38]. Therefore, we can set a constraint on the mean segment length (which is widely used in geology) to a desirable scale for a specific engineering application, e.g., for reservoir engineering, a scale of hundreds of meters will be sufficient. The damage zones distributed at the relevant scale play an essential role in the corresponding fluid flow problem.

We assume the length of a fault trace from a seismic map to be L=10 km. Figure 5 shows different generations of fault segments. They have the mean segment length of 1117 m, 353 m, and 116 m, respectively. For a desired value of the mean segment length, we repeat the simulation until the desired L¯se value is reached within a given tolerance.

Figure 5

Fault segments with different mean segment lengths. Due to self-similarity, (b) is similar to (a) and (c) to (b) after a magnification.

Figure 5

Fault segments with different mean segment lengths. Due to self-similarity, (b) is similar to (a) and (c) to (b) after a magnification.

Real fault segments are not perfectly parallel. We can add randomness in the orientations of the segments and make the geometry closer to reality. Figure 6 shows the tilted fault segments with the mean segment length of 1223 m, 399 m, and 136 m, respectively.

Figure 6

Tilted fault segments with different mean segment lengths. Due to self-similarity, (b) is similar to (a) and (c) to (b) after a magnification.

Figure 6

Tilted fault segments with different mean segment lengths. Due to self-similarity, (b) is similar to (a) and (c) to (b) after a magnification.

Both parallel and tilted fault segments may have crossovers (Figure 7(a)) or intersections (Figure 7(b)) during the subdivision. It is not physically meaningful to have such patterns during the formation of fault segments, and it is also inappropriate to continue the subdivision with crossovers and intersections. Therefore, we implement an algorithm to merge the close and intersecting segments to prevent their crossover and intersections. The merged segments are shown in green in Figure 7. The merge algorithm naturally sets a limit on the mean length of fault segments, because the distance between subfault segments decreases with more subdivisions. Figure 4(b) shows a lognormal fit of the segment length before and after adding the merge procedure to the algorithm. The mean segment length is increased because of the merging. The segment lengths approximately follow the lognormal distribution as expected in Equation (5).

Figure 7

Crossover and intersection of fault segments. Red and blue segments are younger subfaults coming from the their old generations, respectively. Green segments are merged fault segments.

Figure 7

Crossover and intersection of fault segments. Red and blue segments are younger subfaults coming from the their old generations, respectively. Green segments are merged fault segments.

2.2. A Modified Algorithm to Calculate Fractal Dimension

Fractal dimension is a measure of a set’s complexity or its irregularity. Mandelbrot [42] defined a fractal as the following: “a set for which the Hausdorff–Besicovitch dimension strictly exceeds the topological dimension,” where the topological dimension refers to the number of coordinates associated with each element in the set. The example chosen by Mandelbrot [42] deals with measuring the length of the coastline of Britain. In this example, the shorter the yardstick length, r, the longer the length of the coastline Lr is, and L tends to increase without a limit. Mandelbrot [42] suggested that an empirical relationship between L and r is as follows:

According to Mandelbrot [42], D in Equation (6) represents the fractal dimension of the coastline. The fractal dimension characterizes the complexity of the coastline, and a larger D reflects a higher complexity of the coastline.

Mandelbrot [42] also proposed other methods to calculate the fractal dimension of a set. One of the methods is to draw circles of a chosen radius r to cover the set completely using a minimum number N of such circles (to “cover a set” means that the set lies inside of the covering circles. A sketch map is shown in Figure 8(a)).

Figure 8

(a) Illustration of a set composed of one line segment fully covered, partially covered, and not covered by a circle; (b) the circle-based method used to calculate the fractal dimension in Okubo and Aki’s paper [1987]; the blue parts of the line segments are clearly not covered by circles and cause an underestimation of the number of circles Nr.

Figure 8

(a) Illustration of a set composed of one line segment fully covered, partially covered, and not covered by a circle; (b) the circle-based method used to calculate the fractal dimension in Okubo and Aki’s paper [1987]; the blue parts of the line segments are clearly not covered by circles and cause an underestimation of the number of circles Nr.

An example of the coastline of Scotland is shown in Figure 9. The image of Scotland has a size of 509×550 pixels, and the size of 1 pixel is about 1.6 km. In Figure 10(a), we plot the number of covering circles N versus circle radius r. They are plotted as a straight line on a double-logarithmic plot, with the slope, or the fractal dimension, DN=1.24 (R2=0.997). Figure 10(b) shows the dependence of the coastline length, Lr=2Nr, on the circle radius. It has a plateau for large radius values (>500 pixels), which can be related to skipping over the fine details of the coastline. For the radius values below 1 pixel, no finer image details can be resolved. The transitional part between the small and large radii can be fitted by a straight line, whose slope provides the fractal dimension. However, an arbitrary choice of data points for line fitting is subjective and may cause significant variations of the corresponding fractal dimension. Instead, we fit the L vs. r dependence with a scaled logistic function, Equation (7), and then, we find the derivative of this logistic function, Equation (8).
Figure 9

Map of Scotland showing the west coastline in red (original image retrieved from the internet).

Figure 9

Map of Scotland showing the west coastline in red (original image retrieved from the internet).

Figure 10

(a) The number of covering circles versus their radius for the coastline of Scotland. (b) The coastline length versus the covering circle radius for the coastline of Scotland. The red line is the scaled logistic fit of the data. The blue tangent line yields the maximum fractal dimension (these results are from the program written by P. Benito and T. Patzek for a graduate course at UC Berkeley). (c) The derivative of the scaled logistic function.

Figure 10

(a) The number of covering circles versus their radius for the coastline of Scotland. (b) The coastline length versus the covering circle radius for the coastline of Scotland. The red line is the scaled logistic fit of the data. The blue tangent line yields the maximum fractal dimension (these results are from the program written by P. Benito and T. Patzek for a graduate course at UC Berkeley). (c) The derivative of the scaled logistic function.

From the plot of the derivative, Figure 10(c), more accurate fractal dimension is determined as the maximum fractal dimension of the set, denoted as DL. Once we plot the tangent line with a slope of (1DL) (the blue line in Figure 10(b)) and intersect the upper and lower bounds of the total length (green lines in Figure 10(b)), we get the upper and lower limits of the circle radii for the maximum fractal dimension.

Using this methodology, we find that the fractal dimension is 1.48 for the coastline of Scotland. The corresponding upper and lower limits on radius are 6 and 177 pixels, respectively. The two fractal dimensions, DN and DL, are different; however, DL reveals more information about the complexity of a given set and the applicable range of fractal characteristics, while DN averages over all data points.

The circle-based method was adopted by Okubo and Aki [55] to calculate the fractal dimension of the San Andreas fault system. However, from their sketch map (Figure 2 in Okubo and Aki [55] or Figure 8(b)), we see that the circles do not cover their set completely, and the number of circles, Nr, is underestimated.

Therefore, we propose a modified algorithm and calculate the fractal dimension of parallel fault segments. The algorithm is accurate for parallel and vertically well-sorted fault segments. It can be extended to the case of subparallel and not well-sorted fault segments that are common in reality, but with lower accuracy. Figure 11 shows a set of parallel and well-sorted line segments to demonstrate the algorithm. The red circle is the circle used in the previous step to partially cover the line segments. We need to find the next circle and subsequent circles to cover the set with a minimum number of circles. The proposed algorithm minimizes the number of circles using the following two rules:

  • (1)

    A single circle should partially cover as many line segments as possible, i.e., in Figure 11, a single circle should cover all n line segments

  • (2)

    With the same number of line segments covered, the circle should cover as large a length of segment parts which lie inside that circle as possible, i.e., both the blue and green circles cover n line segments. But, we choose the green circle, because it covers a larger length, which is proved in the appendix

Figure 11

A sketch map illustrating the process used to cover a set of line segments completely with a minimum number of circles. The red circle is the circle chosen in the previous step; the next circle that continues to cover the segments with maximum coverage (green or blue) should be determined. The appendix proves that the green circle yields a fuller coverage.

Figure 11

A sketch map illustrating the process used to cover a set of line segments completely with a minimum number of circles. The red circle is the circle chosen in the previous step; the next circle that continues to cover the segments with maximum coverage (green or blue) should be determined. The appendix proves that the green circle yields a fuller coverage.

In Figure 11, if the uppermost and lowermost line segments are covered completely, all other line segments in between will be covered automatically. Therefore, the only situation to be considered is a set composed of two line segments, with as many line segments in between as possible. Intuitively, the blue circle is the first option for covering these line segments, because it passes through the two intersection points, A and B. However, if we rotate the blue circle with respect to point B in a counter-clockwise direction, the coverage of line segment ln, BD, will decrease and the coverage of line segment l1, AC, will increase. If the rotation can make the covered length larger, the maximum length should emerge when point C is exactly on the line segment l1 and BC is the diameter of the circle O2. At this critical point, the covered length is maximum, and beyond it, the covered lengths of both line segments l1 and ln decrease. In the appendix, we demonstrate that the green circle (scenario 2) always covers a larger length than the blue circle (scenario 1). Figure 11 is a more general case. If AB is perpendicular to the line segments, the blue and green circle overlap with each other completely.

With the correct choice of the circle to cover the line segments, the pseudocode of our algorithm is shown in Algorithm 1. The complete program can be found online (doi:10.4121/uuid:38a2ed3e-f456-429c-9fcd-704b392bc8a6).

    Algorithm 1: Covering line segments with minimum number of circles.
  •  begin

  •   Initialization;

  •   while Not all line segments in the set are fully covered do

  •    while Not all line segments in the current circle are fully covered do

  •     Find the next circle centre;

  •     Find the line segments which intersect the circle;

  •     Check whether the line segments are fully covered by current or previous circles;

  •     if All line segments in the current circle are fully covered then

  •      break;

  •     end

  •    end

  •    Remove the fully covered line segments and add new line segments;

  •   end

  • end

In the initialization part of the algorithm, one rotates the set of line segments to make them parallel to the X-axis (please note that all segments are assumed to be parallel or subparallel). Afterwards, the start points (with smaller X-coordinates) of the line segments are sorted along the Y-direction and afterwards along the X-direction. The coverage of the line segments starts from the left-most position of the set and continues rightwards. The method used to find the first circle for each iteration is sketched in Figure 12. For a given radius r, (i) draw a circle with the left-most point, A, in l1 as the center and 2×r as the radius; (ii) find the intersection points on the right side of all intersecting line segments. If there is no other line segment (Figure 12(a)), use AB as the diameter of the first circle. Suppose now that there are more line segments that intersect the green circle, and the intersection points are B,B, and B. Start from the farthest line segment (l3), draw the red circle O, and calculate the covered length CB+AD. If the covered length is larger than or equal to 2×r, circle O is the first circle (Figure 12(b)). Otherwise, find the next farthest line segment and repeat the procedure. This is also the method used to implement rule 1, which makes a single circle partially cover as many line segments as possible.

Figure 12

Two scenarios to find the first circle (red circle) in each iteration.

Figure 12

Two scenarios to find the first circle (red circle) in each iteration.

If the situation shown in Figure 13 happens where EC cannot be covered by the circle O1, we divide the line segment EB into two segments, EC and CB. Segment EC will be a new line segment mentioned in the pseudocode and will be covered in the next iteration. For the last circle in each iteration, circle O2 in Figure 13, the line segments do not span the circle, and we count the number of circles with a fraction, the ratio between the maximum length covered (FN in Figure 13), and the diameter 2r instead of 1. By doing this, we can avoid an overestimation of the number of circles. To make the program more efficient and save computing time, we also divide the domain into blocks, and only neighbouring segments that share the same blocks are checked for coverage in each iteration.

Figure 13

A sketch map of the procedures used to cover the line segments. The coverage will start from the left-most position and continue rightwards. The left part, EC, will become a new segment and be covered in the next iteration. Circle O2 will not be regarded as a whole circle but as a fraction of the circle (FN/2R).

Figure 13

A sketch map of the procedures used to cover the line segments. The coverage will start from the left-most position and continue rightwards. The left part, EC, will become a new segment and be covered in the next iteration. Circle O2 will not be regarded as a whole circle but as a fraction of the circle (FN/2R).

In this section, we implement the modified algorithm to calculate the fractal dimension of natural and generated fault segments.

3.1. Fractal Dimension of Real Fault Maps

The discontinuous and segmented characteristics of faults can be analyzed using fault maps obtained from high-resolution seismic surveys or outcrops. In this paper, three fault maps are chosen to calculate the fractal dimensions using the proposed algorithm. The first one is the San Andreas fault system (Figure 1(b)). The northern part of the San Andreas fault system captures the geometric features that are common to all strike-slip faults with segmentations [22]. The fault segments that are mapped on the surface [50, 56] are oriented NW-SE with the echelon patterns. A sketch map illustrating the process used to calculate the fractal dimension is shown in Figure 14 on the map. The second fault is a Quaternary fault array in the central Apennines in Italy (Figure 2 in [23]). The Monte Vettore fault array (MVFA) in the central Apennines fault system is a part of a Quaternary fault system with the roughly N-S trending left-lateral strike-slip fault zone. The MVFA is characterized by a fault trace pattern dominated by the pervasive linkages among the fault segments of various lengths, making up the geometrically complex, multiscale, minor fault zones within the array. The third fault map is the Top Kharaib fracture lineaments from the Lekhwair Field (Figure 6 in [24]). On top of the Kharaib Formation, there are two well-developed fault families oriented NW-SE and NNW-SSE. These faults are steep, and the fault throws reach a typically visible offset on seismic data with a resolution of approximately 10 meters.

Figure 14

The method of covering the fault from Figure 1(b) by circles with the radius, r=420 m. The red line segments are simplifications of the real fault segments.

Figure 14

The method of covering the fault from Figure 1(b) by circles with the radius, r=420 m. The red line segments are simplifications of the real fault segments.

Considering the San Andreas fault, the fault length is plotted in Figure 15(a) versus the covering circle radii. This figure illustrates that the length of the fault and the circle radius have a linear relationship within a certain range on the log-log plot. After fitting with a scaled logistic function, we find that the fractal dimension is 1.29. The lower limit on the circle radius is about 66 meters, and the upper limit is about 226 meters. Below the lower limit, all details of fault segments are already captured due to the limited resolution of the map. Consequently, the total length of the fault does not change and yields a plateau in the figure. Above the upper limit, all details about the segmentation are ignored, and that makes the total length of the fault remain unchanged. In this case, the fault traces are relatively nonfragmented, simple, and nonfractal geometric sets. For the other two fault maps we chose, the results are shown in Figures 15(b) and 15(c). The fractal dimensions calculated from the fits are 1.29 and 1.34, respectively. The upper and lower limits for fault map 2 are 461 and 2597 meters, respectively. For fault map 3, the corresponding values are about 847 and 3054 meters. A detailed summary of the results for these three real fault maps is listed in Table 1.

Figure 15

(a–c) Fault length versus radius of covering circles. (d–f) The derivative of the scaled logistic function. (a) The fault map 1 (northern part of the San Andreas fault system), (b) the fault map 2 (Quaternary fault array in the central Apennines in Italy), and (c) the fault map 3 (the Top Kharaib fracture lineaments in the Lekhwair Field). The red line is the scaled logistic fit to the data. The blue tangent line yields the maximum fractal dimension.

Figure 15

(a–c) Fault length versus radius of covering circles. (d–f) The derivative of the scaled logistic function. (a) The fault map 1 (northern part of the San Andreas fault system), (b) the fault map 2 (Quaternary fault array in the central Apennines in Italy), and (c) the fault map 3 (the Top Kharaib fracture lineaments in the Lekhwair Field). The red line is the scaled logistic fit to the data. The blue tangent line yields the maximum fractal dimension.

Table 1

Fractal dimensions of three real fault systems.

CaseL¯se (m)Fractal dimensionUpper limit (m)Lower limit (m)
Map 114101.2966226
Map 275171.294612597
Map 394271.348473054
CaseL¯se (m)Fractal dimensionUpper limit (m)Lower limit (m)
Map 114101.2966226
Map 275171.294612597
Map 394271.348473054

3.2. Fractal Dimension of Generated Fault Segments

In this section, the fractal dimensions of the fault segments generated with the algorithm described in Section 2.1 are calculated. We choose to generate fault segments with two segments in each generation, and no rotation of each segment is implemented because the fractal calculation is accurate for parallel and vertically well-sorted fault segments. The generated segments with L¯se=110, 311, and 1664 meters are shown in Figures 16(a)–16(c), respectively.

Figure 16

Fault segments with different mean segment lengths. Due to self-similarity, (b) is similar to (a) and (c) to (b) after a magnification. There are two segments in each generation.

Figure 16

Fault segments with different mean segment lengths. Due to self-similarity, (b) is similar to (a) and (c) to (b) after a magnification. There are two segments in each generation.

The calculated fractal dimension is shown in Figure 17 and summarized in Table 2. The fault length is a linear function of the circle radii within a certain range on a log-log plot. When the radius is small or large, the length of the fault is constant. The lower limit and upper limit of the radius yield an interval in which the fault segments show the fractal or hierarchical self-similar characteristic. The values of these limits are directly related to the resolution of the sample set. The upper and lower limits of the radius increase with an increase in the mean segment length. If we assume that there are two segments with equal length generated in each generation and the step length is k times of the segment length, we can have an exact fractal geometry with the fractal dimension:
Figure 17

(a–c) Fault length versus radius of covering circles. (d–f) The derivative of the scaled logistic function. (a) The generated fault segment 1, L¯se=110 meters; (b) the generated fault segment 2, L¯se=311 meters; (c) the generated fault segment 3, L¯se=1664 meters. The red line is the scaled logistic fit of the data. The blue tangent line yields the maximum fractal dimension.

Figure 17

(a–c) Fault length versus radius of covering circles. (d–f) The derivative of the scaled logistic function. (a) The generated fault segment 1, L¯se=110 meters; (b) the generated fault segment 2, L¯se=311 meters; (c) the generated fault segment 3, L¯se=1664 meters. The red line is the scaled logistic fit of the data. The blue tangent line yields the maximum fractal dimension.

Table 2

Fractal dimensions of three generated fault segments.

CaseL¯se (m)Fractal dimensionUpper limit (m)Lower limit (m)
11101.262.710.5
23111.3111.640.7
316641.2541.5124.4
CaseL¯se (m)Fractal dimensionUpper limit (m)Lower limit (m)
11101.262.710.5
23111.3111.640.7
316641.2541.5124.4

However, the nonuniform step length and width in the generated fault segment can change the fractal dimension, as well as the nonuniform segment length. Therefore, the fractal dimensions of the generated fault segments vary between 1.25 and 1.31. The fractal dimension is slightly smaller than the dimension of real fault segments, because generated fault segments are less complex in geometries than the real fault segments. If more complexities, such as more segments in each generation, orientation rotations are incorporated in the generated fault segments, and their fractal dimension can increase and be closer to reality. In general, the generated fault segments are comparable with natural fault segments from a fractal dimension perspective.

Current research on fault damage zones mostly focuses on the spatial distribution of fault damage zones at outcrops but rarely incorporates fault damage zones in subsurface flow problems, possibly because subsurface fault segments and their associated damage zones are undetectable. Although the method proposed in this paper cannot provide the unique solution of detailed fault segments, it offers possible scenarios, which are consistent with outcrop observations and share similar fractal characteristics. Furthermore, a reasonable structure of fault segments makes the modelling of associated damage zones around each fault segment possible, which is essential for subsurface fluid flow and will be discussed in future research.

In this paper, we have proposed a new method to generate fault segments from an imprecise fault trace obtained from a seismic survey or outcrop. The generated fault segments can be used to analyze the fault damage zone and the fault-related flow problems. These segments have characteristics similar to those found in the natural fault segments observed from outcrops or seismic surveys. The generated segment lengths follow approximately a lognormal distribution. The segment geometry recovers the hierarchical self-similar architecture, which is consistent with the natural fault segments. We also propose a modified algorithm to calculate the fractal dimension of the fault segments. Our algorithm can cover the target set more accurately, and it improves on Okubo and Aki’s method [1987]. The calculated fractal dimensions of the generated and real fault segments are similar. This finding validates our approach from a fractal dimension perspective.

Appendix

Derivation of the Largest Covering Circle

Two cover scenarios are shown in Figures 18 and 19. The distance=h between the line segments l1 and l2 is known, the radius of the circle is R, and the length of O1M is unknown and denoted by a. Here, we prove that for a given h, no matter how a changes, scenario 2 will yield a longer covered length.

Figure 18

A sketch map of scenario 1, where for a fixed red circle, the next circle is blue, and it passes through the two intersection points, A and B.

Figure 18

A sketch map of scenario 1, where for a fixed red circle, the next circle is blue, and it passes through the two intersection points, A and B.

Figure 19

A sketch map of scenario 2, where for a fixed red circle, the next circle is green, and it is obtained by rotating the blue circle with respect to the intersection point A until point N is on the segment l2 (AN is the circle diameter).

Figure 19

A sketch map of scenario 2, where for a fixed red circle, the next circle is green, and it is obtained by rotating the blue circle with respect to the intersection point A until point N is on the segment l2 (AN is the circle diameter).

In Figure 18, the covered length is L1=AX+BY. It is simple to prove that ΔANO2 is equal to ΔBMO1, because the following equations hold
Therefore, NO2=MO1=a, and L1 is
In Figure 19, ADN=π/2 because AN is the diameter of circle O2. ΔO1CF and ΔACD are similar, because the following equations (Equation (A.3)) hold:
Thus, we have
Since
we also have
Since ED=h/tanAEC, we have
Denote AEC=α and O1EC=β. Then,
From ΔO1CF and ΔAED, we also know that
By denoting O2AN=θ, we calculate the covered length in scenario 2:
The difference between L1 and L2 is denoted as f.
Take the derivative of f with respect to a.
dL1/da is trivial.

It is easy to see then that if a<h/2, X1>0.

dL2/da is nontrivial.
The first term is simple, because θ is independent of a.
Based on the chain rule, we calculate
The left hand term in Equation (A.16) is
Combined with Equation (A.9), we have the following:
Therefore, the left term in Equation (A.16) is as follows:

It is obvious that X2 in Equation (A.19) is positive.

The right term in Equation (A.16) is
The left term in Equation (A.20) is
The right term in Equation (A.20) is
Therefore, Equation (A.20) becomes

If a<h/2, X3 in Equation (A.23) is positive.

Therefore, the derivative of fa, Eq. (A.15), is

fa is a monotonically increasing function for ah/2. When a=h/2, f has the maximum value, which is 0. If a>h/2, it is easy to see that X1<0, X2>0, and X3<0. Therefore, fa is a monotonically decreasing function with respect to a>h/2, and the maximum value is still 0 when a=h/2. When ah/2, f<0, which means that scenario 2 always covers more segment length than scenario 1.

The open-source program of calculating the fractal dimension of parallel or subparallel fault segments is available online (doi:10.4121/uuid:38a2ed3e-f456-429c-9fcd-704b392bc8a6).

Key Points. We propose a new method for the automatic generation of fault segments from an imprecise fault trace. We improve on the method used by [55] to calculate the fractal dimension of parallel or subparallel fault segments. We verify that the generated fault segments are similar to natural fault segments in terms of their fractal and geometrical properties.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

This project was supported by the baseline research funding from KAUST to Prof. Tad W. Patzek. The authors would like to thank all editors and anonymous reviewers for their comments and suggestions.

1.
de Joussineau
G.
Aydin
A.
Segmentation along strike-slip faults revisited
Pure and Applied Geophysics
 , 
2009
, vol. 
166
 
10-11
(pg. 
1575
-
1594
)
[PubMed]
2.
de Joussineau
G.
Mutlu
O.
Aydin
A.
Pollard
D. D.
Characterization of strike-slip fault-splay relationships in sandstone
Journal of Structural Geology
 , 
2007
, vol. 
29
 
11
(pg. 
1831
-
1842
)
[PubMed]
3.
Segall
P.
Pollard
D.
Mechanics of discontinuous faults
Journal of Geophysical Research: Solid Earth
 , 
1980
, vol. 
85
 
B8
(pg. 
4337
-
4350
)
4.
De Joussineau
G.
Aydin
A.
The evolution of the damage zone with fault growth in sandstone and its multiscale characteristics
Journal of Geophysical Research: Solid Earth
 , 
2007
, vol. 
112
 
B12
[PubMed]
5.
Flodin
E. A.
Aydin
A.
Evolution of a strike-slip fault network, Valley of Fire State Park, southern Nevada
Geological Society of America Bulletin
 , 
2004
, vol. 
116
 
1
(pg. 
42
-
59
)
[PubMed]
6.
Myers
R.
Aydin
A.
The evolution of faults formed by shearing across joint zones in sandstone
Journal of Structural Geology
 , 
2004
, vol. 
26
 
5
(pg. 
947
-
966
)
[PubMed]
7.
Cartwright
J. A.
Trudgill
B. D.
Mansfield
C. S.
Fault growth by segment linkage: an explanation for scatter in maximum displacement and trace length data from the Canyonlands Grabens of SE Utah
Journal of Structural Geology
 , 
1995
, vol. 
17
 
9
(pg. 
1319
-
1326
)
[PubMed]
8.
Knott
S. D.
Beach
A.
Brockbank
P. J.
Brown
J. L.
McCallum
J. E.
Welbon
A. I.
Spatial and mechanical controls on normal fault populations
Journal of Structural Geology
 , 
1996
, vol. 
18
 
2-3
(pg. 
359
-
372
)
[PubMed]
9.
Aydin
A.
Fractures, faults, and hydrocarbon entrapment, migration and flow
Marine and Petroleum Geology
 , 
2000
, vol. 
17
 
7
(pg. 
797
-
814
)
[PubMed]
10.
Caine
J. S.
Evans
J. P.
Forster
C. B.
Fault zone architecture and permeability structure
Geology
 , 
1996
, vol. 
24
 
11
(pg. 
1025
-
1028
)
11.
Jourde
H.
Flodin
E. A.
Aydin
A.
Durlofsky
L. J.
Wen
X.-H.
Computing permeability of fault zones in eolian sandstone from outcrop measurements
AAPG Bulletin
 , 
2002
, vol. 
86
 
7
(pg. 
1187
-
1200
)
12.
Odling
N.
Harris
S.
Knipe
R.
Permeability scaling properties of fault damage zones in siliclastic rocks
Journal of Structural Geology
 , 
2004
, vol. 
26
 
9
(pg. 
1727
-
1747
)
[PubMed]
13.
Aydin
A.
Berryman
J. G.
Analysis of the growth of strike-slip faults using effective medium theory
Journal of Structural Geology
 , 
2010
, vol. 
32
 
11
(pg. 
1629
-
1642
)
[PubMed]
14.
Childs
C.
Watterson
J.
Walsh
J.
Fault overlap zones within developing normal fault systems
Journal of the Geological Society
 , 
1995
, vol. 
152
 
3
(pg. 
535
-
549
)
[PubMed]
15.
Huggins
P.
Watterson
J.
Walsh
J.
Childs
C.
Relay zone geometry and displacement transfer between normal faults recorded in coal-mine plans
Journal of Structural Geology
 , 
1995
, vol. 
17
 
12
(pg. 
1741
-
1755
)
[PubMed]
16.
Marchal
D.
Guiraud
M.
Rives
T.
Geometric and morphologic evolution of normal fault planes and traces from 2D to 4D data
Journal of Structural Geology
 , 
2003
, vol. 
25
 
1
(pg. 
135
-
158
)
[PubMed]
17.
Nicol
A.
Watterson
J.
Walsh
J.
Childs
C.
The shapes, major axis orientations and displacement patterns of fault surfaces
Journal of Structural Geology
 , 
1996
, vol. 
18
 
2-3
(pg. 
235
-
248
)
[PubMed]
18.
Soliva
R.
Benedicto
A.
Schultz
R.
Maerten
L.
Micarelli
L.
Displacement and interaction of normal fault segments branched at depth: implications for fault growth and potential earthquake rupture size
Journal of Structural Geology
 , 
2008
, vol. 
30
 
10
(pg. 
1288
-
1299
)
[PubMed]
19.
Vermilye
J. M.
Scholz
C. H.
Fault propagation and segmentation: insight from the microstructural examination of a small fault
Journal of Structural Geology
 , 
1999
, vol. 
21
 
11
(pg. 
1623
-
1636
)
[PubMed]
20.
Walsh
J.
Bailey
W.
Childs
C.
Nicol
A.
Bonson
C.
Formation of segmented normal faults: a 3-D perspective
Journal of Structural Geology
 , 
2003
, vol. 
25
 
8
(pg. 
1251
-
1262
)
[PubMed]
21.
Willemse
E. J.
Segmented normal faults: correspondence between three-dimensional mechanical models and field data
Journal of Geophysical Research: Solid Earth
 , 
1997
, vol. 
102
 
B1
(pg. 
675
-
692
)
22.
Aydin
A.
Schultz
R. A.
Effect of mechanical interaction on the development of strike-slip faults with echelon patterns
Journal of Structural Geology
 , 
1990
, vol. 
12
 
1
(pg. 
123
-
129
)
[PubMed]
23.
Cello
G.
Fractal analysis of a Quaternary fault array in the central Apennines, Italy
Journal of Structural Geology
 , 
1997
, vol. 
19
 
7
(pg. 
945
-
953
)
[PubMed]
24.
Filbrandt
J. B.
Al-Dhahab
S.
Al-Habsy
A.
Harris
K.
Keating
J.
Al-Mahruqi
S.
Ozkaya
S. I.
Richard
P. D.
Robertson
T.
Kinematic interpretation and structural evolution of North Oman, Block 6, since the Late Cretaceous and implications for timing of hydrocarbon migration into Cretaceous reservoirs
GeoArabia
 , 
2006
, vol. 
11
 
1
(pg. 
97
-
140
)
25.
Koledoye
A. B.
Aydin
A.
May
E.
Three-dimensional visualization of normal fault segmentation and its implication for fault growth
The Leading Edge
 , 
2000
, vol. 
19
 
7
(pg. 
692
-
701
)
[PubMed]
26.
Benedicto
A.
Rives
T.
Soliva
R.
The 3D fault segmentation development-a conceptual model. Implications of fault sealing
First eage international conference on fault and top seals-what do we know and where do we go?
2003
Montpellier, France
pg. 
97
 
27.
Kattenhorn
S. A.
Pollard
D. D.
Integrating 3-D seismic data, field analogs, and mechanical models in the analysis of segmented normal faults in the Wytch Farm oil field, southern England, United Kingdom
AAPG Bulletin
 , 
2001
, vol. 
85
 
7
(pg. 
1183
-
1210
)
28.
Koledoye
B. A.
Aydin
A.
May
E.
A new process-based methodology for analysis of shale smear along normal faults in the Niger Delta
AAPG Bulletin
 , 
2003
, vol. 
87
 
3
(pg. 
445
-
463
)
[PubMed]
29.
Soliva
R.
Benedicto
A.
A linkage criterion for segmented normal faults
Journal of Structural Geology
 , 
2004
, vol. 
26
 
12
(pg. 
2251
-
2267
)
[PubMed]
30.
Crider
J. G.
Pollard
D. D.
Fault linkage: three-dimensional mechanical interaction between echelon normal faults
Journal of Geophysical Research: Solid Earth
 , 
1998
, vol. 
103
 
B10
(pg. 
24373
-
24391
)
31.
Mansfield
C.
Cartwright
J.
Fault growth by linkage: observations and implications from analogue models
Journal of Structural Geology
 , 
2001
, vol. 
23
 
5
(pg. 
745
-
763
)
[PubMed]
32.
Young
M. J.
Gawthorpe
R. L.
Hardy
S.
Growth and linkage of a segmented normal fault zone; the Late Jurassic Murchison-Statfjord North Fault, northern North Sea
Journal of Structural Geology
 , 
2001
, vol. 
23
 
12
(pg. 
1933
-
1952
)
[PubMed]
33.
Otsuki
K.
Dilov
T.
Evolution of hierarchical self-similar geometry of experimental fault zones: implications for seismic nucleation and earthquake size
Journal of Geophysical Research: Solid Earth
 , 
2005
, vol. 
110
 
B3
[PubMed]
34.
Manighetti
I.
Caulet
C.
De Barros
L.
Perrin
C.
Cappa
F.
Gaudemer
Y.
Generic along-strike segmentation of A far normal faults, East Africa: implications on fault growth and stress heterogeneity on seismogenic fault planes
Geochemistry, Geophysics, Geosystems
 , 
2015
, vol. 
16
 
2
(pg. 
443
-
467
)
[PubMed]
35.
Ando
R.
Tada
T.
Yamashita
T.
Dynamic evolution of a fault system through interactions between fault segments
Journal of Geophysical Research: Solid Earth
 , 
2004
, vol. 
109
 
B5
[PubMed]
36.
Lunn
R. J.
Willson
J. P.
Shipton
Z. K.
Moir
H.
Simulating brittle fault growth from linkage of preexisting structures
Journal of Geophysical Research: Solid Earth
 , 
2008
, vol. 
113
 
B7
[PubMed]
37.
Kim
Y.-S.
Peacock
D.
Sanderson
D. J.
Mesoscale strike-slip faults and damage zones at Marsalforn, Gozo Island, Malta
Journal of Structural Geology
 , 
2003
, vol. 
25
 
5
(pg. 
793
-
812
)
[PubMed]
38.
Kim
Y.-S.
Peacock
D. C.
Sanderson
D. J.
Fault damage zones
Journal of Structural Geology
 , 
2004
, vol. 
26
 
3
(pg. 
503
-
517
)
[PubMed]
39.
Fachri
M.
Tveranger
J.
Braathen
A.
Schueller
S.
Sensitivity of fluid flow to deformation-band damage zone heterogeneities: a study using fault facies and truncated Gaussian simulation
Journal of Structural Geology
 , 
2013
, vol. 
52
 (pg. 
60
-
79
)
[PubMed]
40.
Kolyukhin
D.
Schueller
S.
Espedal
M. S.
Fossen
H.
Deformation band populations in fault damage zone—impact on fluid flow
Computational Geosciences
 , 
2010
, vol. 
14
 
2
(pg. 
231
-
248
)
[PubMed]
41.
Qu
D.
Tveranger
J.
Fachri
M.
Influence of deformation-band fault damage zone on reservoir performance
Interpretation
 , 
2017
, vol. 
5
 
4
(pg. 
SP41
-
SP56
)
[PubMed]
42.
Mandelbrot
B. B.
Fractals: form, chance, and dimension, vol. 706
 , 
1977
San Francisco, USA
WH Freeman
43.
Aviles
C.
Scholz
C.
Boatwright
J.
Fractal analysis applied to characteristic segments of the San Andreas fault
Journal of Geophysical Research: Solid Earth
 , 
1987
, vol. 
92
 
B1
(pg. 
331
-
344
)
[PubMed]
44.
Berkowitz
B.
Hadad
A.
Fractal and multifractal measures of natural and synthetic fracture networks
Journal of Geophysical Research: Solid Earth
 , 
1997
, vol. 
102
 
B6
(pg. 
12205
-
12218
)
45.
Cowie
P. A.
Sornette
D.
Vanneste
C.
Multifractal scaling properties of a growing fault population
Geophysical Journal International
 , 
1995
, vol. 
122
 
2
(pg. 
457
-
469
)
[PubMed]
46.
Hirata
T.
Fractal dimension of fault systems in Japan: fractal structure in rock fracture geometry at various scales
Fractals in geophysics
 , 
1989
Springer
(pg. 
157
-
170
)
47.
Roy
A.
Perfect
E.
Dunne
W. M.
McKay
L. D.
Fractal characterization of fracture networks: an improved box-counting technique
Journal of Geophysical Research: Solid Earth
 , 
2007
, vol. 
112
 
B12
[PubMed]
48.
Thingbaijam
K. K. S.
Nath
S. K.
Yadav
A.
Raj
A.
Walling
M. Y.
Mohanty
W. K.
Recent seismicity in Northeast India and its adjoining region
Journal of Seismology
 , 
2008
, vol. 
12
 
1
(pg. 
107
-
123
)
[PubMed]
49.
Sharp
R.
Implications of surficial strike-slip fault patterns for simplification and widening with depth
U.S. Geological Survey Open-File Report
 , 
1979
, vol. 
79
 
1239
(pg. 
66
-
78
)
50.
Wallace
R. E.
Surface fracture patterns along the San Andreas fault
Conference on tectonic problems of the San Andreas fault system
1973
Stanford, CA, USA
(pg. 
248
-
250
)
51.
Aydin
A.
Page
B. M.
Diverse Pliocene-Quaternary tectonics in a transform environment, San Francisco Bay region, California
Geological Society of America Bulletin
 , 
1984
, vol. 
95
 
11
(pg. 
1303
-
1317
)
52.
Davy
P.
On the frequency-length distribution of the San Andreas fault system
Journal of Geophysical Research: Solid Earth
 , 
1993
, vol. 
98
 
B7
(pg. 
12141
-
12151
)
53.
Bonnet
E.
Bour
O.
Odling
N. E.
Davy
P.
Main
I.
Cowie
P.
Berkowitz
B.
Scaling of fracture systems in geological media
Reviews of Geophysics
 , 
2001
, vol. 
39
 
3
(pg. 
347
-
383
)
[PubMed]
54.
Makarov
P.
Evolutionary nature of structure formation in lithospheric material: universal principle for fractality of solids
Russian Geology and Geophysics
 , 
2007
, vol. 
48
 
7
(pg. 
558
-
574
)
[PubMed]
55.
Okubo
P. G.
Aki
K.
Fractal geometry in the San Andreas fault system
Journal of Geophysical Research: Solid Earth
 , 
1987
, vol. 
92
 
B1
(pg. 
345
-
355
)
[PubMed]
56.
Brown
R. D.
Jr.
Wolfe
E. W.
Map showing recently active breaks along the San Andreas fault between Pt. Delgada and Bolinas Bay, California
1970
Tech. Rep., US Geological Survey
Exclusive Licensee GeoScienceWorld. Distributed under a Creative Commons Attribution License (CC BY 4.0).