## Abstract

Fracture mechanics modelling of fracture pattern development was used to analyse pattern geometry and population statistics for natural opening-mode fractures. Orthogonal fracture network geometries were generated under biaxial extension loading conditions from a slightly anisotropic initial strain state. Fracture statistics were analysed by grouping all fracture orientations into one population for these unique orthogonal pattern geometries. Fracture aperture distributions resembled negative exponential curve shapes, consistent with published observations for stratabound fractures in sedimentary rock. Fracture length distributions had a strongly power-law shape, and showed that longer fractures grew first and reached their fullest extent before shorter fractures began propagating. The power-law shape of the length distribution was first established by the growth of the longest fractures in the population, followed by the later growth of shorter fractures that extended the power-law shape to smaller sizes. The shortest fracture length at which the power-law distribution was truncated varied with the magnitude of the applied strain. Other variations in fracture pattern results were tied to mechanical layer thickness and subcritical crack growth propagation properties of the fractured media.

The characterization of natural, opening-mode fractures is a challenging task, particularly in the subsurface, where only limited information is available. Discovering the important processes responsible for fracture pattern development and the fundamental systematics of fracture attribute populations is essential for fracture analysis. Important fracture attributes include fracture orientation, planarity, spacing, aperture (opening displacement), and length. It is well established from mechanics that opening-mode fractures form perpendicular to the local least compressive principal stress (Lawn & Wilshaw 1975; Pollard & Aydin 1988). For conditions favouring vertical fractures (such as when the maximum compressive stress is vertical), this implies that an isolated opening-mode fracture (free from any mechanical interaction with other nearby fractures) propagates perpendicular to the remote least compressive horizontal stress, σ_{hmin}. Based on both empirical and theoretical analysis, average fracture spacing is often noted to be controlled by mechanical layer thickness (e.g. Ladeira & Price 1981; Narr & Suppe 1991; Bai & Pollard 2000), with the added distinction that the amount of clustering in the pattern may vary as a function of mechanical crack interaction and the degree of pattern development (Delaney *et al.* 1986; Rives *et al.* 1992; Olson 1993; Wu & Pollard 1995; Gillespie *et al.* 2001; Olson 2004). The degree of planarity of a fracture (how straight it propagates) has been interpreted to be primarily controlled by the stress state and the influence of other nearby fractures (Pollard *et al.* 1982; Olson & Pollard 1989; Bai & Gross 1999). For nominally horizontal, bedding plane outcrops, it can be said that crack path curving is enhanced by close crack-to-crack proximity and by a low in-plane differential stress (σ_{hmax}−σ_{hmin}). Straighter crack paths result when fractures are far apart and when the in-plane differential stress is high.

Fracture size attributes such as aperture and length are often analysed with respect to their population statistics. Many researchers have proposed that the cumulative frequency of fracture aperture populations can be described as power-law functions (e.g. Sanderson *et al.* 1994; Clark *et al.* 1995; Marrett *et al.* 1999; Ortega & Marrett 2000). Gillespie *et al.* (1999) found a distinction between stratabound (fracture height contained by mechanical layer thickness) and non-stratabound opening-mode fractures, observing power-law aperture distributions for non-stratabound fractures and non-power-law distributions (normal, log-normal or negative exponential) for stratabound fractures. Ortega & Marrett (2000) also observed an influence of mechanical layer thickness on aperture distributions, but interpreted all distributions as power-law, with different fitting exponents depending on whether fracture length was less than or greater than mechanical layer thickness. Gillespie *et al.* (2001) also observed length distributions to be influenced by fracture height containment, observing lognormal distributions for stratabound and power-law distributions for non-stratabound fractures. Ortega & Marrett (2000) interpreted their fracture length distributions as power-law, but emphasized that length determination is not straightforward for natural fractures, particularly when they can be made up of en echelon segments that may or may not appear linked depending on the scale of observation. Although Ackermann *et al.* (2001) studied experimentally generated normal faults, not opening mode fractures, it is worth noting that they found the nature of fault length data to change depending on the magnitude of the strain, with power-law length distributions at low strain and negative exponential length distributions at higher strain. They also found mechanical layer thickness to be important in controlling fault length and displacement development.

Most of the work on attributes discussed above was limited to singling out a particular set of nominally parallel fractures, even when multiple sets were present, because multiple fracture orientations are typically interpreted to be of different age and thus caused by different deformation events (Pollard & Aydin 1988). However, there are special cases where multiple fracture orientations can result from a single deformational event. Polygonal fracture patterns, such as those found in mudcracks and columnar jointing of volcanic rock, have nominally random orientations of fractures that can be linked to a single causative event, usually an isotropic tension caused by desiccation or cooling (Lachenbruch 1962; DeGraff & Aydin 1993). Orthogonal patterns may also be generated from a single deformational event, and are very commonly found in outcrops of thin sedimentary units. The focus of this study is patterns that are typically referred to as cross-jointed or ladder patterns (Pollard & Aydin 1988; Rawnsley *et al.* 1992; Gross 1993; Rives *et al.* 1994). An important attribute of cross-jointed fracture patterns is that there is typically an older, through-going fracture set that is very straight and regularly spaced (Fig. 1; Gross 1993; Bai *et al.* 2002). Later cross-joints typically propagate perpendicular to the through-going set and abut into them, although there can be a wide variation in the planarity and orientation of these later fractures.

As described by Bai *et al.* (2002), there are multiple mechanisms that can be responsible for the generation of orthogonal fracture sets. The mechanism of interest in this study has been termed stress release, where the tensile stresses caused by a remotely applied extensional strain can be completely relieved by crack propagation perpendicular to that strain (Olson 1993; Wu & Pollard 1995). This stress release builds over time as the fracture pattern develops, and when the initial fracture set becomes well-developed (through-going and regularly spaced), the local principal stress directions can change, causing later fractures to grow in a different (typically orthogonal) direction. Bai & Pollard (2000) illustrated how the reduction in tensile stress in one direction, induced by fracture propagation, can locally cause the principal stress directions to flip by 90°, resulting in a subsequent flip in fracture propagation direction. They identified values of stress anisotropy and the ratio of fracture spacing to mechanical layer thickness for which cross-jointing is expected to occur.

In addition to the stress analysis work exemplified by the study by Bai & Pollard (2000), several studies have been performed where fracture propagation was explicitly modelled. Olson (1997) simulated cross-joint development using an isotropic initial state followed by monotonic uniaxial extension, generating an early set of nominally parallel, through-going fractures between which later cross-joints developed, but the earlier fractures varied considerably in their orientation and were not as universally through-going as often observed in the field. Similar simulations run with biaxial extension loading from an isotropic initial state (Tuckwell *et al.* 2003; Rijken 2005; Olson *et al.* 2007) resulted in more polygonal fracture geometries. Olson *et al.* (2007) also showed that an initial straight fracture propagation followed by cross-jointing can be achieved under biaxial extension given a small initial strain anisotropy, but if that initial anisotropy is too large, it can prevent the principal stress direction change induced by stress release (Bai *et al.* 2002), suppressing the formation of cross-joints.

This study further analyses the systematics of fracture patterns that grow under biaxial extension, conditions that are conducive to orthogonal pattern development. A process-oriented approach is adopted, using mechanics-based numerical experiments to examine the development of pattern geometry and fracture attribute distributions through time. The effects of mechanical layer thickness and subcritical index (stress corrosion cracking) are investigated.

## The fracture mechanics model

Natural, opening-mode fractures have been postulated to initiate from rock imperfections such as microcracks, grain boundary defects, eccentrically shaped pores, fossil fragments and flute casts (Engelder 1987; Pollard & Aydin 1988; Engelder & Lacazette 1990). These mechanical imperfections or flaws can range in size from tens of microns to a few centimeters. The numerical modelling presented here mimicked this initial flawed state of rock by randomly locating vertical, randomly oriented flaws of 0.2 m in length, all centred on a given horizontal plane within a semi-infinite elastic body. An additional restriction on the location of the initial flaws was that their centers could not be closer than two flaw lengths, preventing initial flaw overlap or intersection. Although the 0.2 m flaw length is an order of magnitude or more larger than suggested by natural observations, smaller initial flaws were considered computationally impractical, given that the boundary element size of the numerical simulation scales with the initial flaws (each flaw is sub-divided into two 0.1 m crack elements, to always have two independent crack tips). Fracture propagation was modelled by adding additional boundary elements of constant size (0.1 m length) at the crack tip, so the longer the fracture, the greater the number of elements it contained. For the 20 m by 20 m horizontal dimensions of the modelled fractured region, it would take a fracture with 200 elements of 0.1 m length to completely span the region from side to side. Reducing the initial flaw size by an order of magnitude (1 cm) would result in an order of magnitude more elements to populate fractures of a given length (e.g. 2000 elements for a 20 m fracture). Given the practical code limitation of about 10 000 single precision elements for a computer with 2 Gigabytes of memory, using the smaller flaw or element size would restrict the model to the growth of very few longer fractures (only five or 10). As the goal of the simulations was to generate realistic outcrop-scale fracture patterns containing hundreds of individual fractures, the larger initial flaw size was considered essential. In addition, previous simulation studies (Olson 1993; 2004) have demonstrated that fractures are self-organizing after propagation is initiated, such that the final fracture patterns generated are largely independent of the initial flaw distributions and are more dependent on mechanical boundary conditions and material properties. However, the main point of concern remains that the initial flaw length will influence the strain at which fracture propagation is initiated (the stress at the crack tip scales with the square root of fracture length).

The numerical fracture model integrates the effects of mechanical crack interaction, fracture containment within prescribed mechanical layers (3D effects), subcritical propagation velocities, and displacement-driven loading conditions in a linear elastic media. The fractured region (Fig. 2) is an imaginary layer within a semi-infinite, rectangular parallelepiped (vertical dimension infinite). The fractures are prescribed to fully extend across the imaginary layer from the top to bottom but cannot propagate vertically across the layer boundaries. Thus, the fractures are all centred on the same horizontal plane with respect to *z*, have a constant height equal to the mechanical layer thickness, and can only propagate laterally, as described by Olson (2004). The fractures are assumed to be sharp-tipped at top and bottom (as if the layer boundaries were welded), and thus the size of their stress perturbation scales with their shortest dimension (typically the height for layer-contained or strata-bound fractures) (Olson 2003). This 3D aspect of the mechanics problem is accounted for by incorporating a 3D correction factor (Olson 2004) in the 2D, plane-strain, displacement discontinuity solution (Crouch 1976). This geometric conceptualization is chosen to correspond to the propagation of natural fractures within thin sedimentary units, where the fracture height is constrained by bedding planes between differing lithologies, and the fracture length is typically much greater than the fracture height. A mechanical effect that is not treated in this model is the influence of the mechanical properties of the layers adjacent to the fractured layer (Helgeson & Aydin 1991; Fischer *et al.* 1995), which can be critical in determining whether fractures will be contained within a particular mechanical layer or not. However, once fracture containment is presumed, the layering heterogeneity is considered to be a second-order effect on the mechanics of stratabound fractures (Bai *et al.* 2000).

The model is implemented, as mentioned above, by subdividing fractures into straight, equal-length boundary elements of constant displacement discontinuity (opening and/or shearing) (Crouch 1976). The normal and shearing displacement discontinuities (*D*_{n} and *D*_{s}, respectively) of the numerical elements are referenced to the *x–y* plane and represent the opening (aperture) and shear offset of a fracture at a particular location along its length. Crack propagation direction and velocity depend upon the opening- and shearing-mode stress intensity factors, *K*_{I} and *K*_{II}, which can be computed as a function of the displacement discontinuities at the crack tip element as given by (Olson 1991)
1
where *E* is Young's modulus, ν is Poisson's ratio, and Δ*a* is the crack element length. The primary mechanism for propagation is stress corrosion cracking (subcritical crack growth), where propagation occurs at a stress intensity factor below fracture toughness (*K*_{Ic}) but above some minimum threshold stress intensity factor (*K**_{I} ≈ *K*_{Ic}/10; Segall 1984), typically related to the hydrolytic weakening of intergranular bonds caused by pore waters (Anderson & Grew 1977; Atkinson & Meredith 1987). Given subcritical crack growth conditions (*K**_{I}≤*K*_{I}<*K*_{Ic}), propagation velocity, *V*, can be written as (Atkinson 1984; Swanson 1984; Olson 1993)
2
where *n* is the subcritical index and *A* is a proportionality constant. *K*_{tip} is based on maximum circumferential stress theory (Erdogan & Sih 1963) and can be written as
3
where θ is the counter-clockwise positive angle defining propagation direction, found by solving
4
Straight (in-plane) crack propagation occurs when *K*_{I}>*K**_{I} and *K*_{II}=0, which corresponds to a propagation angle of θ=0° and *K*_{tip}=*K*_{I}. The maximum kink in crack path direction occurs when *K*_{I}=0 and *K*_{II}≠0, resulting in an angle of θ=±70.5° for positive or negative *K*_{II}, respectively.

The displacement discontinuities for a given boundary element, *D*_{n} and *D*_{s}, are calculated as a function of the local stress boundary conditions on the element (Fig. 3a), the remote stress boundary conditions resolved onto the element (Fig. 3b), and stresses caused by the mechanical interaction with other nearby elements (Crouch 1976). The local and remote stress boundary conditions can be summarized with the driving stress term from fracture mechanics (Pollard & Segall 1987), which can be defined for mode I opening as
5

where σ_{n}^{local} is the locally applied normal traction acting on the element (typically an internal fluid pressure in nature) and σ_{n}^{remote} is the resolution of the remote boundary stresses acting perpendicular to the element. For mode II (in-plane) shearing, the driving stress can be defined as
6
where σ_{s}^{local} is the locally applied shear traction acting on the element (typically a frictional resistance when the crack surfaces are in contact) and σ_{s}^{remote} is the resolution of the remote boundary stresses acting parallel to the element. Mechanical interaction with other fracture elements adds additional resolved normal and shear stresses (Crouch 1976).

For the modelling presented here, initial conditions were divided into stress and displacement components. The initial stress conditions were set to give both opening- and shearing-mode driving stresses equal to zero. The condition of Δσ_{II}=0 was accomplished by having an isotropic remote horizontal stress state (σ_{xx}=σ_{yy}), which resolved no shear stress on any vertical fracture element regardless of its strike in the horizontal plane, comparable with the non-tectonic reference state of McGarr (1988). The condition of Δσ_{I}=0 was achieved with an internal fluid pressure in the fracture elements that just balanced the remote horizontal compressive stress (Segall 1984; Olson 1993). Presuming a mechanism of natural hydraulic fracturing, Engelder & Lacazette (1990) proposed that, under static conditions, the fluid pressure in the fracture should be equal to the pore pressure in the formation (*P*_{p}). Consequently, these initial stress conditions could be used to represent any depth in the subsurface where the pore pressure was equal to the horizontal *in situ* stress (*P*_{p}=σ_{xx}^{remote}=σ_{yy}^{remote}). Prescribing the exact loading or burial path leading to this initial state is beyond the scope of this paper, but there are various models in the literature that could be manipulated to reach this point (Prats 1981; Zoback & Healy 1984; Blanton & Olson 1999). In addition, a zero shear-stress condition was imposed on the vertical boundaries (perimeter) of the solid to allow it to freely expand or contract in response to the lateral normal displacement loading (making these boundaries principal stress planes).

The displacement boundary conditions were applied to the perimeter of the fractured body (Fig. 2) to simulate the effects of horizontal tectonic strains. An initial shortening strain in the *x*-direction of ε_{xx}=−3.6×10^{−5} (ε_{yy}=0) was assumed to be the cumulative effect of horizontal tectonics up to the beginning of the simulations, during which no crack growth had occurred other than the creation of the starter flaws. The value of the initial ε_{xx} was chosen to be large enough to align initial crack propagation in the *x*-direction (Olson *et al.* 2007), avoiding random polygonal fracture pattern development such as modelled by Tuckwell *et al.* (2003), but small enough to prevent the purely straight crack propagation indicative of large horizontal stress anisotropy (Olson & Pollard 1989; Olson *et al.* 2007). Fracture propagation was induced by the application of additional horizontal tectonic displacements constituting a biaxial extension of ε_{xx}=ε_{yy}=9.1×10^{−5} applied at a strain rate of 1.2×10^{−19} s^{−1} over a time period of *c*. 24 Ma (Fig. 4). The total extension was small in magnitude but is comparable with that indicated by the observation of joints in the field (Segall 1984). The strain rate was chosen to represent a relatively stable sedimentary basin, not an active tectonic margin. Based on geodetic and geological data, crustal strain rates have been estimated to range from 1×10^{−13} to 1×10^{−20} s^{−1} (Kreemer & Holt 2001; Kreemer *et al.* 2002). The horizontal extension that drove crack propagation was achieved in the numerical model by imposing 10 discrete increments of horizontal normal displacement at regular time intervals to the vertical edges of the semi-infinite solid, representing a constant strain rate deformation in a step-wise fashion. The tectonic (remote) strain anisotropy was quantified with the strain ratio, ε_{D}, defined as
7

ε_{D} >1 indicates loading that is more anisotropic than uniaxial conditions (ε_{D}=1.0), whereas ε_{D}<1 indicates loading that is less anisotropic than uniaxial conditions. The limit of ε_{D}=0.0 represents isotropic conditions. Thus, the larger the value of ε_{D}, the stronger the horizontal tectonic strain anisotropy, and the stronger the elastic-induced, tectonic stress anisotropy to align cracks for straight propagation (Olson & Pollard 1989).

## Numerical modelling results

Simulations were performed with 200 initial flaws, mechanical layer thicknesses of 2 and 4 m, and subcritical index values of *n*=20, 40 and 80. The other mechanical properties used in the simulations are indicated in Table 1. The subcritical index values were based on the experimental range of values typical for sandstone, which average around *n*=50 (Holder *et al.* 2001; Rijken *et al.* 2002). The horizontal dimensions of the semi-infinite body were 22 m×22 m, and the body perimeter was divided into 2 m long boundary elements. Fracture propagation was limited to a 20 m×20 m sub-region centred within that body to avoid unwanted numerical edge effects caused by fracture interaction with these large edge elements (see Olson 1993).

### Orthogonal trace pattern geometry

All simulations had the identical strain history (Fig. 4): the strain state was strongly anisotropic (ε_{D}>1) until about 9 Ma, while ε_{xx} was contractional; it was uniaxial (ε_{D}=1) for the duration of one loading step, while ε_{xx}=0 (from 9 to 11 Ma); and it became weakly anisotropic (ε_{D}<1) after 11 Ma, when ε_{xx} became extensional (positive) and no longer hindered cross-fracture growth. The results (Fig. 5) showed a progression in fracture pattern geometry from a clearly cross-jointed pattern, with the through-going fractures aligned with the initial *x*-directed shortening strain (Fig. 5a), to a more polygonal pattern with a slight preferential fracture trend in the *x*-direction (Fig. 5f). Increasing layer thickness and increasing subcritical index had similar effects on the fracture pattern: they both tended to increase fracture spacing and increase the tendency of crack paths to curve. The increase in fracture spacing with mechanical layer thickness was a direct result of the increasing size of the stress shadow around the fracture (area of stress release), which suppressed later fracture growth in the vicinity of early formed fractures and scaled in size with the fracture height (equal to mechanical layer thickness) (Pollard & Segall 1987; Bai & Pollard 2000; Olson 2004). This effect was clearly demonstrated by the difference between the patterns with a 2 m layer thickness (Fig. 5a, c and e) and those with a 4 m layer thickness (Fig. 5b, d and f).

Fracture paths were less straight for the thicker mechanical layer (*t*=4 m) than for the thinner one (*t*=2 m), all other factors being equal. This was caused by the higher crack-induced stresses in the thicker layer case that increased the degree of mechanical interaction between cracks. The magnitude of crack-induced stresses for a mode I crack should scale roughly with the magnitude of the stress intensity factor, *K*_{I}, at the tip (Lawn & Wilshaw 1975). This scaling is precise in the near-tip region (at distances from the crack tip that are much less than the fracture dimensions) and only approximate further away from the crack, but it is sufficiently accurate to allow better understanding of the change in crack interaction with changing mechanical layer thickness. A strata-bound fracture can be approximated as an elliptical fracture with a height (minor axis) equal to *t* (mechanical layer thickness) and a length (major axis) of *L*. The solutions for the minimum and maximum stress intensity factors around an elliptical fracture are (Kanninen & Popelar 1985)
8
acting at the ends of the major axis (the lateral tips), and
9
acting at the ends of the minor axis (the top and bottom of the fracture). *Ei*(*k*) is the complete elliptic integral of the second kind, and *k*^{2}=[(*L*/2)^{2}− (*t*/2)^{2}]/(*L*/2)^{2}. Based on the numerators of equations (8) and (9) (the variation of *Ei*(*k*) with *t* is weak) and assuming a given length and constant driving stress, both (*K*_{I})_{min} and (*K*_{I})_{max} increase strongly with increasing mechanical layer thickness, *t*. Given that the remote differential stress (caused in these simulations by the remote strain anisotropy) will not vary with layer thickness, the overall effective mechanical interaction between cracks is enhanced by the higher crack-induced stresses in thicker mechanical layers, and this mechanical interaction leads to more crack path curving.

The effect of increasing subcritical index on the fracture pattern spacing and fracture path curving was attributed to fracture propagation timing. Fracture growth and strain anisotropy v. time for the 2 m layer cases (Fig. 6) showed that fracture propagation (indicated by the increasing number of crack elements) began later for higher subcritical index cases. This delay in crack propagation with increasing subcritical index was an outgrowth of subcritical propagation velocity law (equation (1)). Because (*K*_{I}/*K*_{IC}) <1 for subcritical crack growth, the higher the subcritical exponent, *n*, the lower the propagation velocity for cracks of a given stress intensity factor. In the numerical model, crack growth occurred only when enough time had elapsed to add an entire boundary element (0.1 m long in this case). Consequently, there was a characteristic time delay proportional to the propagation velocity before crack length increased. For instance, the *n*=20 case had no recordable growth until about 6.2 Ma into the simulation. For *n*=40, it took 11.5 Ma for growth to start, and for *n*=80 it took 17.8 Ma.

The increased delay before fracture propagation began for the higher subcritical index material affected fracture spacing in a couple of ways. First, at any given point in the strain history, the flaw population had a range of propagation velocities that was a function of stress intensity factor according to equation (1). Only the fastest growing cracks accumulated enough length to add crack elements. Given the finite time span of the simulations, the lower the subcritical index of the material, the more the total fracture length accumulated (see the number of elements added for various subcritical indices shown in Fig. 6), and consequently the lower the average fracture spacing. The minimum average spacing between the through-going, *x*-parallel fractures was about 2.7 m for a subcritical index of *n*=20 (Fig. 5a). Increasing the subcritical index to *n*=40 increased average fracture spacing to about 3.1 ms (Fig. 5c), and at *n*=80, average fracture spacing increased to 4.4 m (Fig. 5e). Similar fracture spacing dependence on subcritical index was obtained for a layer thickness of 4 m (Fig. 5b, d and f), with values of 3.7, 4.4, and 7.3 m for subcritical indices of *n*=20, *n*=40, and *n*=80.

Another subcritical index related factor that influenced the spacing between the more through-going, *x*-parallel fractures was the point in the strain history at which cross-fractures started forming. Once cross-fractures began to form, it was more difficult for *x*-parallel fractures to increase in length without intersecting other fractures. For the patterns generated in this study, cross-jointing always occurred after the total *x*-strain became extensional (positive), or in other words, when ε_{D}<1. For the 2 m thick layer with *n*=20, *x*-parallel propagation of through-going fractures occurred between 6.2 and 10.7 Ma (Fig. 7a and b), with 1502 crack elements added while ε_{D}≥1.00 (Fig. 6). Cross-fracturing began at 11.3 Ma (Fig. 7c), with predominantly diagonal trends to the pre-existing, *x*-parallel fractures. At 13 Ma (Fig. 7d), cross-fracturing continued, but predominantly orthogonal to the through-going fractures. Some subsequent in-filling to the cross-fractures also appeared in the form of younger and shorter *x*-parallel fractures, which continued between 16.2 and 22.8 Ma (Fig. 7e). Additional *y*-parallel cross-fractures and some diagonally oriented fractures occurred at this late time as well.

For the 2 m thick layer with *n*=40, fracture propagation did not initiate until about 11.5 Ma (Fig. 6), under remote conditions that were already more isotropic than uniaxial strain (ε_{D}<1). This reduced anisotropy made it easier for cross-fractures to form as a result of the stress relief caused by the early *x*-parallel fracturing. The *x*-parallel fractures were still fairly straight, but the spacing was broader when cross-jointing began than for the *n*=20 case (one less through-going fracture propagated). The crack growth reconstruction for the *n*=40 case (not shown) indicated that cross-jointing began at 19.7 Ma (ε_{D}=0.50), which, compared with the *n*=20 case, was later in the strain history but earlier in the fracture propagation history (only 1231 elements of straight crack growth had occurred prior to cross-jointing with *n*=40, compared with 1502 for *n*=20). The difference between the amount of parallel fracture infilling prior to cross-jointing for *n*=20 v. *n*=40 was even more significant for the case of a 4 m thick layer (Fig. 5b v. Fig. 5d).

For *n*=80 with the 2 m thick layer, propagation initiation was delayed until almost 17 Ma. The initial fracture propagation followed fairly non-planar paths, indicating strong mechanical interaction between nearby fractures (Olson & Pollard 1989). Mechanical interaction was stronger for this case of delayed fracture initiation for two reasons: (1) the strain anisotropy had become weak (ε_{D}=0.57) by 17 Ma, so the remote differential stress that could have aligned cracks into straight paths was weak when cracks started growing; (2) the magnitude of overall strain at the time of crack initiation was higher than for the lower subcritical index cases (Fig. 4), meaning larger stress intensity factors when cracks started growing, and larger crack-induced stresses. The *n*=80 case for the 4 m thick layer shows similar results but with even less *x*-parallel fracture set development.

### Fracture aperture development

Based on field data, Gillespie *et al.* (1999) proposed that aperture distributions for stratabound veins followed normal, log-normal or negative exponential distributions, whereas non-stratabound vein apertures followed power-law distributions. Data from the *n*=20 case for a 2 m thick layer were plotted with best-fit curves for various distribution shapes (Fig. 8). The data did not show power-law character but were more similar to the negative exponential shape, consistent with the findings of Gillespie *et al.* (1999) for stratabound fractures. Looking at the aperture data from all of the simulations (Fig. 9), many exhibited similarities to a negative exponential shape.

In addition to the distribution shape, there were several systematic variations in aperture worth noting (Fig. 9). The 4 m thick layer cases had peak apertures almost twice those from the 2 m thick cases, demonstrating the aperture dependence on layer thickness in the simulations. There was also a systematic variation with subcritical index, where the *n*=80 cases had the greatest number of larger aperture fractures and the fewest smaller aperture fractures (the shallowest slope in the semi-log plot and the lowest ‘*y*-intercept’), whereas the *n*=20 cases had the fewest larger aperture fractures and the most smaller aperture fractures. Thus, as exemplified in the tracemap results (Fig. 5), fewer fractures grow in the higher subcritical index cases, but those fractures that do grow have larger apertures because they have to accommodate the same amount of overall strain and they have fewer nearby fractures that would interfere with their opening via mechanical interaction. (These generalizations exclude the fracture apertures that are less than 0.001 m, which plot as an almost vertical line along the cumulative frequency axis. They represent the flaws that did not grow, and were not considered in the analysis of distributions as they do not reflect the self-organizing growth and mechanical interaction process.)

### Fracture length analysis

Whereas the aperture distribution shapes varied depending on subcritical index values, the length distributions for these same simulations were all strongly power-law (Fig. 10). Regardless of layer thickness or subcritical index, the longer end of the length distributions (greater than 1–2 m in length) plotted as a straight line on a log–log plot with a power-law exponent of *c*. −2. The main difference between length populations for the different cases was that the higher the subcritical index and the greater the layer thickness, the larger the portion of the initial flaw population that did not grow. For the 4 m thick layer, only 15% of the fractures exceeded the initial flaw length (0.2 m) for the *n*=80 case, whereas the *n*=20 case had almost 50% of the fractures growing beyond 0.2 m. Decreasing the layer thickness from 4 to 2 m for *n*=20 increased the percentage of fractures growing beyond 0.2 m from 50% to 80%.

To understand how the power-law length distributions developed, it is interesting to look at the distribution shape through time, using the example of the 2 m thick layer with *n*=20 (Fig. 11). As shown in Figure 6, fracture growth did not start until about 6.2 Ma, so the length distribution at 2.2 Ma simply reflects the flaw population. There was a small amount of growth indicated by 8.4 Ma, including the appearance of the longest fracture in the simulation (open square plotted at 20 m), which completely spanned the breadth of the fractured region in the *x*-direction. Although other early fractures appear to span the entire fractured body (Fig. 7a), most are made up of overlapping en echelon segments, similar to fractures that have been documented in the field (Segall & Pollard 1983; Vermilye & Scholz 1995). The sequence of population distributions (Fig. 11) shows the growth of successively shorter fractures through time, filling in the power-law trend from longest to shortest. This sequence suggests that once a fracture started to grow, it quickly grew to the fullest extent possible, upon which its propagation was arrested. The maximum length of younger fractures was dictated by the probability of intersection with earlier propagated fractures. As time progressed and the fracture pattern filled in, the maximum likely length for a newly propagating fracture diminished. The shorter fractures, thus, were not those whose growth was terminated early in the deformation history, they were the ones that started growing the latest.

## Discussion

There have been several recent papers describing how cross-jointed patterns might develop (e.g. Gross 1993; Rives *et al.* 1994; Bai *et al.* 2002); this study adds the perspective of controlled numerical experiments in fracture generation, along the lines of the study by Tuckwell *et al.* (2003) on polygonal fracture patterns. The modelled patterns (Fig. 5) demonstrated the potential influence of rock with different subcritical indices and varying mechanical layer thickness, indicating that cross-jointed patterns are more likely to form in thinner-bedded units and for low to intermediate subcritical index (*n*<80). The results also showed how a small initial strain anisotropy can guide early fracture propagation direction to establish the through-going set, whereas initial isotropy gives a more polygonal pattern (Tuckwell *et al.* 2003; Olson *et al.* 2007), and a stronger anisotropy suppresses cross-jointing (Bai *et al.* 2002; Olson *et al.* 2003). Consequently, there is a fairly narrow range of boundary conditions from which cross-jointed patterns will evolve, as delineated by Bai *et al.* (2002) based on stress analysis.

The distribution analysis of the numerical fracture patterns presented here is somewhat unusual compared with what is typically done with field data. The analysis in this study lumped fractures of multiple orientations together into one population. Typically, analysed outcrops either have only one fracture orientation or the population is segregated into parallel sets. For instance, the numerical pattern in Figure 5a (as well as the natural patterns in Fig. 1) could be argued to consist of an early, regularly spaced, systematic set of parallel fractures, and a non-systematic assemblage (the cross-joints and later generations) that is probably much younger. Based on this interpretation of systematic v. non-systematic, the cross-joints would justifably be excluded from detailed analysis because they would be presumed to be caused by a different event from that which caused the systematic set. However, knowing *a priori* that the simulated fractures in this study were all generated as a result of one deformational event, it seems clear that all orientations should be analysed as members of one population. The corollary is that similar natural fracture sets should be treated likewise. The scaling of the numerical results, particularly for the lengths, validate this approach in that systematic results were obtained. The development of the fracture pattern through time as displayed in Figure 7 depicts a deformation process of sequential infilling (Gross 1993) in response to monotonically increasing strain, not a series of separate fracture events of varying orientation.

The length distributions generated with the studied loading conditions were clearly power-law in character, but this power-law result does not indicate a universal characteristic for all fracturing processes. Simulations performed with the same numerical model under strongly anisotropic conditions have produced single orientation fracture sets with negative exponential length distributions (Olson 2004). The reason for the difference in distribution type for the cross-jointed v. parallel fracture patterns may be a difference in propagation arrest mechanisms. For the cross-jointed sets (natural and numerical), virtually every fracture abuts into a pre-existing one, suggesting that length is determined by how likely a fracture is to intersect pre-existing ones, and those intersections are typically at right angles. For a population of subparallel, straight fractures, pattern development takes on a different character, with propagation arrest caused by the indirect growth interference of mechanical interaction between non-intersecting, overlapping en echelon fracture tips (Pollard *et al.* 1982; Olson & Pollard 1989). A key characteristic for straight fracture growth, given a non-intersecting mechanism for propagation arrest, is that the application of additional strain can result in the lengthening of previously arrested fractures. This was demonstrated in growth experiments for parallel normal faults by Ackermann *et al.* (2001), where additional strain increased the lengths of all active faults, including the longest ones. For the cross-jointed patterns of this study, additional strain does not lengthen pre-existing fractures but initiates new ones. Thus, the numerical results suggest there is not a universal population scaling law for opening mode fracture lengths, but the scaling varies depending on the nature of the fracture interactions that limit length.

The length development sequence of Figure 11 also brings up an interesting question with regard to how population distributions change as strain accumulates. The power-law length distributions for each snapshot in time were truncated at the shorter end, but became less truncated with further pattern development caused by increasing applied strain. In field data, the change in slope of a power-law distribution toward smaller values is usually attributed to a sampling limitation termed truncation bias (Ortega & Marrett 2000). If the numerical models are a valid proxy for natural fracture pattern development (at least in cross-jointed patterns), truncation bias in natural fracture populations could have a genetic cause related to the amount of strain experienced at a given locality. This is analogous to the concept applied to fracture spacing analysis by Rives *et al.* (1992) and Wu & Pollard (1995), among others, where fracture spacing distributions are interpreted to change through time as a function of how well-developed or ‘saturated’ a particular pattern becomes.

Finally, projecting the power-law distribution of Figure 11 back toward a cumulative frequency of one can be used to estimate the expected abundance of shorter fracture lengths in the pattern if the loading were to continue past 25 Ma. The predicted value of the shortest fracture is about 1.3 m, assuming all flaws eventually propagated and the power-law length distribution was maintained. Clearly, field observations do not support an absence of fractures below 1.3 m, given the common observations of microfractures down to the grain scale (Kranz 1983; Laubach 1997). The documented progression of length distribution through time suggests that part of the solution to this discrepancy is to increase the number of initial flaws; that is, as fracture pattern development involves the infilling of the pattern with the propagation of new fractures, there need to be sufficient initiation sites or flaws for that progression to continue. In defence of the flaw number used for this study (200), the outcrop patterns depicted in Figure 1 do not suggest that further fracture propagation would be required beyond what is exemplified in Figure 5.

Another question that needs to be addressed with respect to flaws is their initial lengths. To increase the initial number of flaws per unit area for a simulation, flaw length would need to be reduced to keep flaws from interacting strongly or amalgamating. The use of shorter, more realistic sized flaws will have a tendency to postpone fracture propagation because of the reduced stress intensity factor for shorter fractures given the same loading. Using equations (8) and (9) and assuming a penny-shaped crack, with *t*=*L*=*r* (where *r* is the radius of the penny-shape), the stress intensity factor can be shown to scale as *K*_{I}α. The results from increasing the subcritical index to *n*=80 (Fig. 5e and f) indicated that delaying propagation initiation under the prescribed loading scenario had the effect of increasing crack path curving and increasing fracture spacing (reducing the amount of propagation). Therefore, if the results of this study were to be applied to a situation where the starter flaws were smaller, it is expected that to obtain a cross-jointed pattern with very straight through-going fractures, a larger initial strain anisotropy would be required to suppress crack path curving and to postpone cross-joint initiation until the through-going fractures became more closely spaced.

## Conclusions

Multiple fracture orientations can be a consequence of multiple loading cycles with different principal stress–strain directions, such as could be caused by different tectonic events. If the earlier fracture sets are sealed with cement, later fractures may propagate across them rather than abut into them, as described by Gillespie *et al.* (2001) for Carboniferous limestones of The Burren, Ireland, where younger joints cross-cut earlier calcite-filled veins. This paper examines something different: the generation of multiple fracture orientations (in the form of a cross-jointed pattern) caused by a single loading cycle. The loading could be a result of uplift and erosion, differential expansion and contraction, or a small strain tectonic event. However, making the connection between the loading character (boundary conditions) and a particular type of geological event was not the purpose of this paper. The goal here was to associate a particular set of mechanical boundary conditions with a fracture pattern. Although this association is not necessarily unique (i.e. there may be multiple sequences of loading that could give similar-looking fracture patterns), and the investigation here was certainly not exhaustive, it does provide some interpretative guidelines for what type of mechanical boundary conditions or loading history are implied by fracture patterns observed in the field. This study showed how layer thickness and subcritical crack index exerted fundamental controls on fracture spacing, fracture timing, and the planarity of fracture paths. Fracture pattern geometry was also strongly influenced by strain anisotropy, which changed in magnitude with time given the conditions of a small initial anisotropy followed by biaxial extension. The trends observed from the simulations presented here show strong similarities to observed field relations, despite the simplifying assumptions required to make numerical modelling feasible, indicating that further investigations, particularly where some limiting assumptions could be relaxed, could provide additional advancement in the understanding of fracture processes.

## Acknowledgments

This study was supported by a grant from the Chemical Sciences, Geosciences and Biosciences Division, Office of Basic Energy Sciences, Office of Science, US Department of Energy program ‘Predicting fracture porosity evolution in sandstone’, and by the industrial associates of the Fracture Research and Application Consortium.

- © The Geological Society of London 2007