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 [1–3]. 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, 4–6]. 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, 9–12].

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 [14–21]. The discontinuous and segmented fault architecture has also been described by high-resolution seismic surveys [22–28].

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 [14–16, 29]. Linkage occurs when the segments overlap during a substantial fault interaction through their associated stress fields [30–32]. 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 [39–41]. 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, 43–48]. 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.

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⁠)

(1)L¯se=2.16S0.89R2=0.87
  • (2)

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

(2)L¯st=0.30S0.77R2=0.84
  • (3)

    Step length x and step width y

(3)x=2.69y0.97R2=0.93

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).
(4)x¯=0.30S0.77,L¯se=2.16S0.89,nse×L¯se−nse−1x¯=L,xi=2.7yi,

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/nse∑i=1nseli (⁠li is the length of the subfault i⁠); the mean step length, x¯=1/nse−1∑i=1nse−1xi (⁠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.

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).
(5)Ln=L×X1×X2×⋯×Xn,

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].

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.

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.

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).

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:
(6)Lr=Ar1−D.

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)).

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).
(7)fx=a+b−a1+edx+c,(8)f′x=−b−adedx+c1+edx+c2.

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 (⁠1−DL⁠) (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

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.

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.

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.

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.

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.

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:
(9)D=ln2ln2−k.

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.

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
(A.1)∠ANO2=∠BMO1=π2,AO2=BO1=R,∠NAO2=∠MBO1.
Therefore, NO2=MO1=a⁠, and L1 is
(A.2)L1=2R2−a2+R2−h−a2.
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:
(A.3)∠ACD=∠ACD,∠O1FC=∠ADC=π2.
Thus, we have
(A.4)O1FAD=O1CAC⟶O1FAD=O1CO1C+R⟶ah=O1CO1C+R⟶O1C=Rah−a.
Since
(A.5)∠AO1E=∠ACE+∠O1EC,∠O1EA=π−∠AO1E2,∠AEC=π−∠O1AE−∠ACE,
we also have
(A.6)∠AEC=π2+∠O1EC2−∠ACE2.
Since ED=h/tan∠AEC⁠, we have
(A.7)ED=x=htan∠ACE2−∠O1EC2.
Denote ∠AEC=α and ∠O1EC=β⁠. Then,
(A.8)ED=x=htanα2−β2.
From ΔO1CF and ΔAED⁠, we also know that
(A.9)cosα=O1FOC=ah−aRa=h−aR,sinβ=aR.
By denoting ∠O2AN=θ⁠, we calculate the covered length in scenario 2:
(A.10)L2=4Rcosθ−x=4Rcosθ−htanα/2−β/2.
The difference between L1 and L2 is denoted as f⁠.
(A.11)f=fa=L1−L2.
Take the derivative of f with respect to a⁠.
(A.12)dfda=dL1da−dL2da.
dL1/da is trivial.
(A.13)X1=dL1da=−2aR2−a2+2h−aR2−h−a2.

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

dL2/da is nontrivial.
(A.14)dL2da=d4Rcosθda−dhtanα/2−β/2da.
The first term is simple, because θ is independent of a⁠.
(A.15)d4Rcosθda=0.
Based on the chain rule, we calculate
(A.16)dhtanα/2−β/2da=hdtanα/2−β/2dα/2−β/2dα/2−β/2da.
The left hand term in Equation (A.16) is
(A.17)dtanα/2−β/2dα/2−β/2=1+tan2α/2−β/2=1cos2α/2−β/2=2cosα−β+1.
Combined with Equation (A.9), we have the following:
(A.18)cosα−β=h−aR2R2−a2+aR2R2−h−a2.
Therefore, the left term in Equation (A.16) is as follows:
(A.19)X2=hdtanα/2−β/2dα/2−β/2=2hh−a/R2R2−a2+a/R2R2−h−a2+1.

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

The right term in Equation (A.16) is
(A.20)dα/2−β/2da=12dαda−12dβda.
The left term in Equation (A.20) is
(A.21)12dαda=12darccosh−a/Rda=−121−h−a/R2−1R=12R2−h−a2.
The right term in Equation (A.20) is
(A.22)12dβda=12darcsina/Rda=121−a/R21R=12R2−a2.
Therefore, Equation (A.20) becomes
(A.23)X3=dα/2−β/2da=121R2−h−a2−1R2−a2.

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

Therefore, the derivative of fa⁠, Eq. (A.15), is
(A.24)dfda=X1+X2×X3>0.

fa is a monotonically increasing function for a≤h/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 a≠h/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.