Quantified fracture (joint) clustering in Archean basement, Wyoming: application of the normalized correlation count method

Fracture spatial arrangement is a key aspect of the structural heterogeneity of the brittle crust (Bear et al. 2012; Laubach et al. 2018a). Understanding the spatial arrangement of opening-mode fractures (joints and veins) has great practical importance. For example, fracture-enhanced permeability is important in aquifers (De Dreuzy et al. 2001; Henriksen & Braathen 2006; Maillot et al. 2016), waste repositories (Barton & Hsieh 1989; Carey et al. 2015), hydrothermal systems (Sanderson et al. 1994; Cox et al. 2001; Hobbs & Ord 2018) and petroleum reservoirs (Laubach 2003; Gale et al. 2014), including those in basement rocks (Cuong & Warren 2009; Belaidi et al. 2018; Bonter et al. 2018). Typical of fracture-controlled permeability is tremendous heterogeneity and anisotropy in fluid flow (e.g. Halihan et al. 2000; Solano et al. 2011). Although fracture network permeability is governed by many factors including aperture and length distributions and connectivity, and by whether fractures are open or sealed, a common cause of permeability variability is heterogeneous spatial arrangement. Here we use a recently introduced method to describe fracture spatial arrangement in a large exposure of texturally uniform Precambrian granitic rock in Wyoming.

‘Spatial arrangement’ describes fracture positions, including whether fractures are clustered, randomly placed or evenly spaced (Laubach et al. 2018a). In conventional usage, spacing is the perpendicular distance between two fractures in the same set (Priest & Hudson 1976; Ladeira & Price 1981; Narr & Suppe 1991; Gillespie et al. 1993), but the prevalence of highly clustered fractures is increasingly being recognized, and what constitutes a useful quantitative measure of clustering is a matter of current interest (Questiaux et al. 2010; Roy et al. 2014; Sanderson & Peacock 2019). For a set of co-planar fractures, a common data collection method is to measure spacing along a straight line of observation (1D scanline) orthogonal to fracture strike. The normalized correlation count (NCC) method (Marrett et al. 2018) accounts for sequences of spacings, including non-neighbour spacings. The technique was implemented with the freely available software CorrCount, accessible via the open-access Marrett et al. (2018) paper. We compare results with conventional spacing analysis methods: descriptive statistics and the ratio of standard deviation to mean or coefficient of variation (Cv).

Here we demonstrate statistically significant self-organized clustering over a length scale range from 10⁻² to 10¹ m in joints in texturally homogeneous Archean granitic rock in a setting distant (>2 km) from folds or faults. These joints are currently conduits for water flow, and their patterns may be useful guides to the spatial arrangement of fractures capable of transmitting fluids in other, less accessible fractured basement rocks. Our example illustrates how NCC can quantify spatial arrangement over a wide range of length scales, measure the pattern and degree of clustering, and test the statistical significance of the results. We use the same method to assess spatial arrangement patterns in six other previously published fractured granite datasets, finding a range of degrees of clustering. Ours is the first application of the NCC technique to fractures in granitic rocks. Results from our example and analysis of datasets from the literature show that opening-mode joints in granites have patterns that range from power-law clustering to indistinguishable from random.

The Teton Range in NW Wyoming is a NNE-trending normal-fault-bounded block of Precambrian crystalline rocks and dominantly west-dipping Cambrian–Mesozoic sedimentary rocks (Love et al. 1992; Roberts & Burbank 1993; Smith & Siegel 2000) (Fig. 1). To the east the range is bounded by the east-dipping Teton normal fault (Love et al. 1992), a structure having Cenozoic–recent movement (Smith et al. 1993; Byrd et al. 1994; Leopold et al. 2007). Within the range are north- to NW-striking reverse faults of probable early Tertiary (Laramide) age (Love et al. 1978; Lageson 1987).

Basement rocks in the Teton Range (Bradley 1956; Love et al. 1992; Reed & Zartman 1973; Zartman & Reed 1998) include complexly deformed Archean gneiss, migmatite and metasedimentary rocks of the Wyoming Province (Chamberlain et al. 2003), Late Archean, mostly unfoliated, granitic intrusive rocks (Mount Owen Quartz Monzonite: Reed & Zartman 1973), and east–west-striking diabase dykes. Our data are from the upper part of Teton Canyon, in granitic rocks mapped as Mount Owen Quartz Monzonite (Reed & Zartman 1973; Love et al. 1992) (Fig. 1). The unit is a silica-rich peraluminous leucogranite (Bagdonas et al. 2016; Frost et al. 2018) forming a discordant pluton (2.55 Ga Mount Owen batholith: Frost et al. 2018), with margins marked by irregular bodies and dykes of pegmatite and aplite, and angular wall-rock inclusions (xenoliths). The exposures we studied are unfoliated, medium to fine, texturally uniform light-coloured granite comprising quartz, microcline and sodic plagioclase with trace biotite and muscovite.

The smoothly undulating, nearly complete exposures we studied retain a fine, striated glacial polish, except when adjacent to some joints where post-glacial erosion has removed the polished surface (Figs 1–5). Locally visible are SE- and NW-facing crescent-shaped chatter mark (fracture) arrays caused by subglacial deformation and linear features on glaciated surfaces caused by post-glacial to recent rock fall (Fig. 2a, b). The glacial and recent geomorphic history of the range is summarized by Foster et al. (2010).

Our fracture orientation and relative timing measurements were made in an outcrop that is c. 0.1 km² (Fig. 1). The scanline was measured with a compass and 30 m plastic tape (centimetre graduations). We orientated the scanline 259° in a direction that maximized the extent of a vegetation- and debris-free outcrop where measurements could be made at approximately uniform elevation (Fig. 3). The scanline is normal to the strike of a main, steeply dipping joint set, and our analysis focuses on these fractures. Scanline location was based on the largest extent of well-exposed rock, not to maximize the number of intersections. Because fractures we measured have long trace lengths relative to outcrop size, moving the scanline around or adding scanlines within the window of good exposure would yield no significant differences in the resulting pattern. Although, across the outcrop, the strike of this set locally has as much as a 10° range, the configuration of the scanline relative to fracture strike means that Terzaghi (1965) corrections of spacing are small. Measurement imprecision of scanline length and fracture placement arises from compliance of the tape, non-uniform tension under field conditions and the non-planar, undulating character of the outcrop (Figs 3 and 4). We judge these errors to be small compared to the scale of the pattern captured in the scanline. We measured kinematic aperture (opening displacement) for representative fractures using a comparator (Ortega et al. 2006) and measured or visually estimated trace lengths for representative fractures.

The NCC technique (Marrett et al. 2018) accounts for the sequence of fracture spacings along a scanline. NCC uses distances between all pairs of fractures, including non-nearest neighbours. The technique provides a quantitative analysis of the degree to which fractures are clustered, and can distinguish between even (or periodic) spacing, clusters arising due to random arrangement and clustering that is stronger than a random signal. NCC is based on the correlation sum, or the two-point correlation function (Bonnet et al. 2001), that calculates the proportion of fracture pairs in a set, including pairs of non-neighbouring fractures, separated by a distance less than each given length scale, λk, in a logarithmically or linearly graduated series of length scales. A correlation count assigned for a given λk is defined as the fraction of all fracture pairs for which the pair's spacing falls between λk + m and λk − m (Marrett et al. 2018), essentially the difference between the correlation sum of λk + m and that of λk − m.

We used the NCC computer program, CorrCount, which provides analytical and Monte Carlo solutions for randomized input spacings and a 95% confidence interval constructed for the randomized sequence (Marrett et al. 2018). The frequencies are normalized against the expected frequency for a randomized sequence of the same fracture spacings at each length scale. Where a length scale's corresponding correlation count falls outside the upper or lower confidence limits, the corresponding fracture spacing can be interpreted to be statistically significant.

We also report conventional fracture spacing statistics (Table 1) and a standard measure of spatial heterogeneity: the coefficient of variation, Cv (e.g. Gillespie et al. 1999). Cv is σ/μ, where σ is the standard deviation of spacings and μ is the mean.

Our outcrop primarily contains opening-mode fractures, including quartz veins (Fig. 5b), apparently non-mineral lined and locally open joints (Figs 4–6), and fractures associated with subglacial chatter marks (Fig. 2). One east–west-striking right-lateral fault with 4 cm displacement was found. The quartz veins have opening displacements (widths) of from <0.1 mm to as much as 1 mm. They are fully sealed with quartz and lack porosity, and have no cement textures evident at hand-specimen scale. These veins have long traces relative to their widths. Aspect ratios locally approach those of the joints, although some veins are markedly wider and shorter than joints. Although locally abundant in the outcrop, veins are rare along the scanline trace. Veins dip steeply and strike NNW, north and ENE. Along the scanline, veins differ in strike from the north-striking joint set we measured; elsewhere in the SW quadrant of the exposure, north-striking joints are subparallel to veins. The quartz veins along the scanline are cross-cut by joints and, in general in this outcrop, veins are the oldest fractures. Veins are not included in our spacing analysis.

Joints are common, although there are large areas of the outcrop that lack fractures of any type (Fig. 3). Joints cross-cut and are readily distinguished from veins. Joints have narrow opening displacement, typically near or below values that reliably can be measured with the Ortega et al. (2006) comparator and hand lens. They are mostly <0.1 mm wide (Fig. 4d). Although joint apertures may have been increased by exhumation or exposure, our measurements put an upper limit on cumulative aperture sizes along the scanline. Many joint traces are marked by dark seams that are likely to be clay-mineral or iron-oxide fills. Some of these have a faintly foliated appearance at the hand-lens scale. These textures may mark minute shear offset. Joints otherwise have no visible mineral fill. Some joints are surrounded by narrow (5–20 cm) irregular halos of red iron oxide stain, marking past fluid flow. In addition, locally along joints, direct evidence of modern fluid flow is apparent (Fig. 5d).

Such narrow fractures are visible mainly because of their great trace lengths. Although most joint traces are long relative to the outcrop size, traces visible at hand-lens scale range from a few to many tens of centimetres to fractures that cross the entire outcrop (>50 m). Short joint traces are undersampled in our 1D analysis. Near our scanline, trace length distributions are censored by outcrop size. About 1 km SE of our outcrop, one ENE-striking fracture zone in basement rocks visible on Google Earth (centred at 43° 41′ 38″ N, 110° 50′ 40″ W) extends for 2 km and can be traced upsection into overlying Paleozoic carbonate rocks. In our outcrop, fracture heights cannot be systematically measured, but fracture traces and outcrop topography imply heights of tens of metres or more for the longest fractures and fracture zones.

All joint sets in our outcrop have steep dips. Although fracture strike dispersion across the outcrop is considerable (Fig. 6), at least four preferred joint strikes are present: NW, north–south, ENE, and WNW to east–west (Fig. 6: sets 1–4). In addition, in nearby exposures, Love et al. (1992) mapped ENE-striking fracture zones (Fig. 6c). The prevalence of various joint orientations varies across our outcrop. Focusing on fractures near our scanline, we designate the two most prominent strike directions set 1 (NW strikes: blue traces in Fig. 5) and set 2 (north strikes: red traces in Fig. 5). For all sets, many fractures are discontinuous and do not intersect other fractures. Where fractures intersect, abutting, stepping and crossing relationships (e.g. Hancock 1985) suggest that east–west-striking joints are locally late. Some abutting relationships, as we note below, imply set 1 fractures predate set 2. Although sets 1 and 2 joints are close to vertical, they have a slight tendency to dip steeply east to ENE (Figs 5a and 6a). Otherwise, the relative timing between the joint sets is mostly ambiguous.

Many joints in both sets 1 and 2 have isolated traces, or are in patterns where two or more fractures have parallel traces over great distances (many metres). En echelon and left- and right-stepping patterns are also apparent. Overlapping tip traces for macroscopic fractures tend to be straight rather than sharply hooked, but close inspection of individual traces shows that they are composed in part of overlapping, gently curved and locally hooked segments, marking small fractures that linked to form longer fractures (e.g. Olson & Pollard 1989; Lamarche et al. 2018). Traces consequently are not everywhere fully interconnected along their intersection with the outcrop surface, making length definition ambiguous, and in some instances the apparent continuity depends on the observation scale. Some en echelon set 1 fractures are connected at their tips by set 2 fractures (Fig. 5e), resulting in a local zigzag pattern of fracture occurrence.

Wing-crack arrays are common and important components of brittle deformation (Willemse & Pollard 1998). Although no offset markers are evident, some set 1 joints are likely to have minute fault displacements based on set 1 in parent configurations relative to set 2 wing cracks (Fig. 5c, d). Fractures in wing-crack configurations (set 2 relative to set 1 in Fig. 5c, d) abut against and extend at a shallow angle from the tips of other fractures (set 1 in Fig. 5c, d). Thus, some set 1 fractures striking between 285° and 310° have north-striking set 2 splays that abruptly diverge from near their tips at angles of 40°–60°. This pattern can arise by right slip on weak planes, such as pre-existing joints. Wing-crack arrays with patterns compatible with left slip are present on some fractures that strike NNE (Fig. 6). Such patterns are likely to mark slip along pre-existing fractures from which the splays diverge (e.g. Willemse & Pollard 1998). The wing cracks propagate in the extensional quadrant of the fault tip. In our pavement, we found little corroborating evidence of slip such as widespread striations (one example was observed) and no macroscopic evidence of shortening, such as stylolites, in the contractional quadrants. Based on the configurations, we interpret set 1 to predate set 2, with some set 2 joints formed as wing cracks, implying a component of right slip on set 1 joints. The concentration of wing cracks and other connecting fractures is apparent at a wide range of scales. Linking set 2 fractures account for locally dense fracture occurrences (Fig. 5e).

Overall, the prevalence of wing-crack arrays, fracture trace patterns that resemble Riedel shear configurations (R1, R2, X) (e.g. Tchalenko 1970; Dresen 1991), smaller fractures near joint traces at acute angles to the main trace, and faintly foliated, pinnate texture and fine-scale undulations along many joint traces, and, locally, striations, suggest that many joints in this outcrop have been sheared. Possibly, populations of pre-existing joints, as well as some quartz veins, were sheared in a consistent pattern compatible with north–south shortening or east–west extension (Fig. 6). Unsheared north-striking opening-mode set 2 joints may have formed as part of this deformation.

Our scanline was aligned to cross set 2 because this set is normal to the long dimension of the outcrop, allowing the longest scanline. Our spatial analysis focused on this set. Set 2 fractures consist of isolated, parallel fractures and fractures in wing-crack arrays. Fracture strain is low (Table 1). We did not assess the heterogeneity of strain (i.e. Putz-Perrier & Sanderson 2008) because for these uniformly very narrow joints we do not have systematic aperture size population data.

Near joints, the glacially polished surface is commonly eroded (Fig. 3), compatible with weathering and fluid flow along open joints. However, joints are also present in glacially polished areas. Here, joint traces appear to abut polished surfaces. We interpret joints to be truncated – or cut by – intact polished surfaces (Fig. 5b–d). If this interpretation is correct, the joints predate glaciation. This inference of relative timing is consistent with the lack of parallelism between joint orientations and current topography, and the lack of evidence of joint concentrations or alignments relative to past ice-flow directions marked by glacial striations. The youngest fractures in the outcrop are chatter marks due to glacial action, readily distinguished from joints and veins (Fig. 2). Also present are linear surface marks (artefacts) caused by post-glacial weathering and erosion (notably, rock fall). These subglacial and late features characteristically damaged the glacial polish that marks much of the surface of the outcrop.

We measured 420 fractures over a scanline distance of 180.4 m. For set 2, spacings have a wide range, from <0.005 m to more than 34 m. Average set 2 spacing for the entire scanline is 0.43 m, but values are strongly skewed towards narrow spacings (Fig. 3c). A measure of spacing heterogeneity is the coefficient of variation (Cv), where Cv is σ/μ: σ is the standard deviation of spacings and μ is the mean. For the overall scanline, Cv is high, 4.86. Occurrence v. distance (‘stick’) plots show anomalously closely spaced fractures separated by large areas where fractures are sparse or absent (Figs 3 and 7). Figure 7 shows clusters at expanded scale, and Figures 8 and 9 show normalized intensity and NCC plots for the scanlines as a whole and for clusters A and B. Descriptive statistics and Cv values for the entire scanline and for subdivisions of the scanline are summarized in Table 1 and in Figures 8 and 9.

CorrCount outputs include normalized fracture intensity (Fig. 8a) and NCC (Fig. 8b) plots. A useful procedure for interpreting CorrCount outputs is to first inspect the normalized intensity plot and then analyse the NCC diagram. On the one hand, for highly clustered fractures, the intensity plot reliably but qualitatively marks regions with fracture concentrations. On the other hand, the value of the NCC plot is such that the degree of clustering can be examined quantitatively across a range of length scales. The scales that can be examined depend, of course, on the dataset; however, for outcrop and subsurface fractures and faults, data ranges are typically from millimetres or less to several hundred metres.

We first inspect the normalized intensity plot. For our scanline, two major clusters within the scanline are apparent (Fig. 8a). Several statistically significant peaks above the upper 95% confidence limit occur at c. 30, 50, 145 and 175 m from the beginning of the scanline. The intensity signals for the rest of the fractures along the scanline lie well within the 95% confidence interval and their distribution is indistinguishable from random. The most statistically significant trough (dips beneath the lower 95% confidence limit) is from 50 to 140 m. Minor troughs occur at the beginning of the scanline, from 0 to 25 m, and around 44, 150 and 165 m.

A stepwise procedure for interpreting NCC diagrams is as follows (Fig. 8b). By comparison with Marrett et al. (2018, fig 12), determine which category of NCC pattern is present. Recognizing that Figure 8b resembles the ‘fractal cluster’ pattern of Marrett et al. (2018, fig. 12h), we note the following features for our dataset (1–5 in Fig. 8b), discussed in order of increasing length scale. Figure 8b, 1: in the NCC plot, a broad, statistically significant, elevated section extends from submetre length scales to c. 21 m, at which point it falls beneath the upper confidence limit and intersects the randomized data curve with a spatial correlation of 1 at near 24 m. The elevated section of the NCC curve indicates that spacings at these length scales are more common than they would be for randomized spacings. Because these are small spacings, clustering is implied. In this section the NCC values decrease as a power law with increasing length scale on the logarithmic NCC plot, suggesting a fractal pattern of spacings within the clusters. The slope and intercept of the best-fit line are characteristics of the pattern.

In the NCC plot (Fig. 8b), the far left of the diagram could be a high-level plateau at small length scales (c. <0.1 m of length scale), indicating that the spacing of the fractures inside the cluster is not organized, although our interpretation is that the pattern marks a power-law decay. The elevated section exhibits multiple peaks averaging around a slight and gradual negative slope before the roll-off point at around 15 m, where the curve rapidly decreases to cross the black line for randomized data. The subsequent trough occurs centred on 54 m. The curve remains under the randomized data curve with a spatial correlation of 1 between 24 and 100 m, the latter exceeds the half-length of the scanline. The range of length scales is small to make a reliable interpretation and close to the point where NCC cannot be calculated because there are no spacings that fit within the window of length scales over which NCC is calculated.

In Figure 8b, 2: the elevated section of the NCC eventually crosses the NCC = 0 value at a length scale of c. 24 m, marking the point at which spacings are comparable with the randomized data (Marrett et al. 2018). This length scale of 24 m is a measure of cluster width. In other words, the NCC intercept near 24 m correlates with the width of the single clusters. In Figure 8b, 3: NCC values for spacings at length scales of >24 m (up to 100 m) lie well below the NCC = 0 line, forming a trough centred on 54 m. This zone marks the gap between clusters where spacings at those length scales are less common than for the randomized data.

The distance between cluster A and cluster B is marked by the right edge of the major trough where the curve intersects the randomized data curve at around 100 m (cf. L in Fig. 8a). Figure 8b, 4: this peak is a measure of cluster spacing, and corresponds to the distance between cluster A and B (Fig. 8a). Because the clusters themselves have a width (rather than a peak at a single length scale), the peak at point 4 will also tend to be quite broad (from 100 to c. 150 m) and in some datasets there may be several harmonics if multiple clusters are observed (cf. Marrett et al. 2018, fig. 12d, f, h).

Although cluster spacing can be detected using logarithmic graduations of length scales, plots of NCC with linear graduations for the entire scanline more readily reveal whether the pattern contains evidence of ‘regularly-spaced fractal clusters’. NCC peaks marked below for logarithmically graduated Figure 9b and d were obtained using linear graduations of length scales because of the scarcity of length scales at high values when using logarithmic graduations. Linear plots of length scale also show more clearly the spacing patterns within the clusters. Figure 8b, 5: nested peaks repeating at increasing length scale are a mark of clusters within clusters. In this example, this signal is noisy, but the fact that there are apparently nested peaks over approximately three orders of magnitude, overall showing a power-law decay in NCC, suggests a fractal pattern. Smaller variations, like those represented by 5 in Figure 8b can be an artefact of the number of length scales chosen for plotting. This is apparent if one increases the number of length scales used to calculate NCC; small variations tend to disappear but the trend above the 95% confidence interval does not.

In order to assess the robustness of the technique for a given dataset it is necessary to inspect the randomness curves: the upper and lower confidence limits. At some length scales they ‘pinch’ in towards the 0 spatial concentration line. This marks the half-length of the scanline (Figs 2 and 8b). NCC values beyond this point are compromised because less than half of the scanline is available for counting the number of fracture spacings at those length scales. The limit is the length of the scanline itself. The overall length of the scanline governs which part of the spacing signal is interpretable. In our dataset, the cluster spacing is wider than the half-length of the scanline because the clusters, A and B, are encountered near the beginning and end of the scanline. The effect of this configuration is to accentuate the amplitude of the trough and peak at 3 and 4, respectively.

Normalized correlation count analysis can be applied to parts, or subdivisions, of the dataset (Fig. 9). In the normalized intensity plot for cluster A (Fig. 9a), four neighbouring, yet distinct, statistically significant peaks occur near 26, 30, 32.5 and 37.5 m, which form a subcluster within cluster A. Two neighbouring peaks occur near 46 and 47.5 m, with the former being the largest of all peaks in the plot, and form the other subcluster within cluster A. Major troughs showing intervals where very few fractures were identified occur along the scanline near the beginning of the scanline, and at 43 and 50–70 m.

The NNC plot for cluster A has a broadly elevated section from submetre length scales to c. 3.5 m, indicating a fractal spatial organization within clusters, as for the whole scanline but with a weaker signal, and at a smaller length scale representing clusters within clusters (nested clusters) discussed for Figure 8b, point 5. The series of peaks and troughs that pass above and below the randomized data confidence limits, from 3.5 to c. 21 m, represent distances between the peaks and subclusters in Figure 9a. The curve dips below NCC = 0 at length scales beyond 21 m because this is longer than the distance between the first and last subcluster (26–47.5 m): that is, the width of cluster A.

In the normalized intensity plot for cluster B (Fig. 9c), peaks are not as distinct and statistically significant as that in cluster A. The highest peak occurs nears 174 m, close to the end of the scanline. A neighbouring peak at around 168 m is also statistically significant. At the 143–147 m and 158–163 m intervals, there are two clusters of peaks that are very narrow, barely above or below the 95% confidence limit, which indicates weak clustering at those intervals. Four troughs occur at near 142, 150, 159 and 171 m. Figure 9d may represent fractal clustering close to indistinguishable from random, although the plot lacks a well-defined low NCC trough (point 4 in Fig. 8b) between cluster width (point 2 in Fig. 8b) and cluster spacing (obtained using linear graduations of length scale). The area above the 95% confidence interval is not continuous, and where it is above the 95% confidence interval it is only for a small range of length scales. The pattern overall is suggestive of fractal clusters but we rate it as nearly indistinguishable from random (Table 1).

Ehlen (2000) described 1D fracture spacing from a large number of different granites. We used data reported by Ehlen for CorrCount analysis. Results of our analysis for a subset of Ehlen's datasets are listed in Table 1 and examples are illustrated in Figure 10.

That the fractures along our scanline are unevenly spaced and occur in clusters is evident from visual inspection (Figs 3 and 4). Anomalously closely spaced fractures in narrow arrays are called fracture swarms or corridors (Laubach et al. 1995; Questiaux et al. 2010; Gabrielsen & Braathen 2014; Miranda et al. 2018; Sanderson & Peacock 2019), but how ‘anomalously’ closely spaced or how narrow or sharply defined an array needs to be to constitute a corridor, and what constitutes a swarm or corridor boundary, have been qualitative assessments. Fracture patterns that are indistinguishable from random will show some degree of clustering (e.g. Laubach et al. 2018a, fig. 1a). The spacing histogram for our dataset (Fig. 3c) merely shows that very closely spaced fractures are common (mean spacing of 0.43 m), an uninformative result that is one of the drawbacks of using frequency distributions to interpret spacing patterns (e.g. McGinnis et al. 2015). NCC provides a way to rigorously identify and quantify corridors.

The ratio of standard deviation of spacings to the mean, the coefficient of variation (Cv), increases with increasing irregularity of spacing, since periodic spacing produces Cv values of zero (Gillespie et al. 1999, 2001) and random locations produce a Cv near 1 (Hooker et al. 2018). For our scanline, the Cv is 4.86, compatible with markedly uneven spacing (Table 1). The marked contrast in spacing within and outside clusters increases the standard deviation. Our Cv value is higher than those reported for fractures in some other granites (Table 1). Ehlen (2000) reported mean spacings and standard deviations for six scanlines in several granites with a mean Cv of 1.41. The Cv for her longest scanline, however, is close to the value we found for the cluster A segment of our scanline. Inspection of the spacing patterns for the Ehlen (2000) scanlines suggests, however, that the Cv, as expected, is not fully documenting the variations in clustering within these examples. The averaged values do not capture mixtures of highly clustered and more regularly space patterns along the scanline.

Our results are the first to use NCC to describe fracture spatial arrangement in granitic rocks. Direct NCC comparison with published NCC results for fractured granites is impossible, but in Table 1 we compare the Wyoming example with NCC values we calculated for spacing datasets in granite reported by Ehlen (2000). Table 1 also reports the Cv, a measure of spacing heterogeneity. Although only recently introduced, the NCC method has been used on fractured carbonate rock outcrops (Marrett et al. 2018), tight gas sandstone horizontal core and image-log data and outcrops (Li et al. 2018), and faults in sedimentary rocks near detachment faults (Laubach et al. 2018b). Compared with other examples, our dataset is more clustered than those in most sedimentary and crystalline rock examples described so far (Table 1).

Interpreting our NCC curve using the Marrett et al. (2018) terminology, we find that of the eight characteristic patterns (their fig. 12), ours resembles that of fractal clusters, where locally close spacings contain smaller self-similar spacing patterns. NCC patterns for the entire scanline (Fig. 8a, b) show a systematic power-law variation of spatial correlation with length scale which suggests that self-organized clustering arose across a wide range of scales. Within the elevated section, repeated patterns of peaks, troughs and plateaus of similar width on the log scale (Fig. 8b, point 5) mark fractal clustering.

Above the 95% confidence line, the slope of the power-law pattern (i.e. exponent of power law) provides a scale-independent measure of the degree of clustering (Marrett et al. 2018, fig. 12), with steeper slopes indicating more intense clustering (i.e. more fractures in clusters at the expense of intercluster regions). Here the slope is relatively high compared to other datasets (Table 1). This means more clustering at small length scales, suggesting there are clusters within clusters, a result supported by the analysis of individual clusters A and B (Fig. 9).

Similar, but weaker, patterns are evident for subdivisions of the scanline within the two main clusters (Fig. 9). Cluster spacing within cluster A is marked by the peaks at about 3 m length scale. Cluster spacing within cluster B is marked by the peaks at about 5 m length scale. Multiples of the dominant spacings are marked by additional peaks separated by intervening troughs (6 m for cluster A; 10 and 15 m for cluster B). Qualitatively, the flatter, more plateau-like within-cluster patterns in cluster B suggest that the organizational style differs. Overall, this subpattern is nearly indistinguishable from random, which may mean that the entire dataset has fractal clustering that is less intense than, for example, in the Pedernales dataset (Marrett et al. 2018).

Variations in the pattern from one cluster to the next might be compatible with self-organized clustering modified locally by some external forcing mechanism, such as the localization of wing-crack arrays by pre-existing set 1 fractures. Clusters also have a range of shapes, from uniform ‘corridor’-like widths to those that vary markedly in width along strike. CorrCount can take fracture size (and other attributes) into account, but we lack the length or aperture measurements that would allow this analysis. Joints here have a narrow range of aperture sizes and small differences in length compared to the scale of the exposure; the spatial organization is likely to reflect the sequence of fracture spacings. The scale of our analysis has limits imposed by finite scanline length and, possibly, limited resolution for small fractures.

The spacing data reported by Ehlen (2000) allow us to calculate NCC values for her crystalline rock examples (Table 1; Fig. 10). Results show a range of patterns. Many of the arrangements are indistinguishable from random, but some have evidence of hierarchical clustering (Fig. 10a, b). Our example and those of Ehlen (2000) show opening-mode joints in granites having patterns that range from indistinguishable from random to power-law clustering. This comparison shows that our Wyoming example fits within the spectrum of clustering found in these other granitic rock scanlines.

Our example and the observations in Ehlen (2000) show why large exposures and methods sensitive to spacing hierarchy are sometimes needed to adequately document the spatial arrangement. The longest scanline in Ehlen (2000) was 74 m and the rest are 40 m long or less. Yet, Ehlen (2000) noted that along her scanlines – qualitatively – patterns varied from regularly spaced to clustered. Because our outcrop is exceptionally long (>180 m) and free of obscuring cover, marked differences in fracture intensity and spacing patterns are apparent. If only part of the scanline were exposed, even in fairly large outcrops, this heterogeneous pattern would be obscured. For example, although from 0 to 50 m fractures are highly clustered and have an NCC pattern that resembles that of the scanline as a whole, from 160 to 175 m the pattern appears more regularly spaced. From 60 to 140 m the rock is largely unfractured; this interval gives a misleading view of the overall intensity. Together with Ehlen's (2000) scanline data and spacing patterns, and evidence in other granite exposures (e.g. Escuder Viruete et al. 2001), both long scanlines and sensitive spatial analysis methods are needed to unravel heterogeneous, hierarchical patterns. The Cv does not account for sequences of spacings, and so is an inadequate measure of spacing irregularity or pattern and cannot capture this hierarchy of patterns.

Correlation count has advantages over other methods that use a cumulative departure from uniform. For example, the correlation sum and Kuiper's method using maximum departures from a uniform distribution (e.g. Sanderson & Peacock 2019) do not reflect the extent to which the pattern is non-random at that location. Although these values can be computed with CorrCount, Marrett et al. (2018) recommended the parameter correlation count because the count allows a signal to be captured at each length scale chosen, rather than up to a given length scale. Correlation count provides more information because it tells the analyst how much more/less common spacings are than random at a given length scale. Thus, cluster spacings and widths, and different patterns at different scales for the same dataset, can be identified. Another advantage of the Marrett et al. (2018) method implemented in CorrCount is the ability to normalize the correlation integral or correlation count to a random rather than a uniform distribution, with the random distribution calculated analytically or modelled with a Monte Carlo method.

Parallel, anomalously closely spaced fractures can arise from several processes that depend on mechanical properties and loading histories (Olson 2004). Without clear evidence of fracture relative timing or spacing variation with strain (i.e. Putz-Perrier & Sanderson 2008; Hooker et al. 2018), the origin of clustering in our outcrop cannot be pinned down. In unstratified rock, stress concentrations near fracture tips can stimulate the development of nearby fractures (Segall & Pollard 1983; Olson 2004). Numerical models show that propagating fracture tips cause nearby offset microfractures to grow, producing tabular, clustered patterns (Olson 1993, 2004). Such a mechanism is compatible with the segmented pattern of fracture traces in our outcrop.

This process of growth and interaction could lead to self-organization. For granitic rock, probably possessing high subcritical crack index values (e.g. Nara et al. 2018) and comprising a thick mechanical unit, densely fractured, widely spaced clusters could be expected. Another process that could produce clusters are simple fault zones formed as oblique dilatant fractures (splay fractures or wing-crack arrays) that link non-coplanar faults side to side and end to end. Closely spaced set 2 fractures connecting tips of set 1 fractures could have formed in this way (Fig. 5d). A similar pattern of progressive joint and joint-fault development has been described from other granitic rocks in the western USA (Martel 1990; Mazurek et al. 2003).

The origins of joints in the Wyoming outcrop are ambiguous. Fractures primarily accommodating opening displacement propagate along a plane of zero shear stress in isotropic rock; specifically, the plane perpendicular to the least compressive principal stress (Pollard & Aydin 1988). Unfortunately, these outcrops provide little basis for determining the timing of joint formation or for inferring their causative stress fields. Cross-cutting relationships with quartz veins show that joints formed well after consolidation of the granite, and overlap by glacial polish is compatible with joints formed prior to glaciation. Steep east to NE set 2 dips could result from joint formation in rocks prior to post-Paleozoic gentle (c. 10°–15°) westwards tilting of the range (Roberts & Burbank 1993; Leopold et al. 2007), although the dip pattern is not well defined and could in any case result from other factors in an igneous mass of unknown dimensions.

These joint sets can be classified as ‘regional’ (Hancock 1985) since they are in otherwise undeformed rocks more than 2 km from the nearest map-scale folds and faults (Fig. 1). Contrasting set orientations and relative timing information imply progressive or episodic joint development with shifting stress orientations. Wing-crack arrays are consistent with some early formed set 1 joints reactivating in right slip. The lack of slip direction indicators and uncertain relative ages of fractures having right and left slip means that the kinematic compatibility of the ensemble of fractures is uncertain.

Without evidence of fracture age, this information still leaves many possibilities for the joint origins. A wide range of loading paths might contribute (e.g. Engelder 1985), including thermoelastic contraction implicated in the formation of some joints in other western USA uplifts (English & Laubach 2017). The shearing of the joint pattern is compatible with north–south shortening or east–west extension. Deformation in these brittle rocks might result from mild bending. The range experienced tilting, and subtle open folds (flexures) are evident elsewhere in the western Tetons. Regional correlates for such deformation are uncertain. Attested regional Late Cretaceous–Eocene shortening directions are east–west, NE (065°) and, later, broadly east–west extension (Powers 1982).

Joints formed in the currently active Cordilleran extensional stress province (Heidbach et al. 2018) and kinematically compatible with the active generally north- and NNE-striking Teton normal fault would be expected to strike north–south to north to NNE. Recent deformation might also be responsible for slip on older NW- and NE-striking joints. On the other hand, joint sets in the granite broadly overlap with the strikes of undated joints and veins in overlying Paleozoic rocks. Temperature patterns of fluid-inclusion assemblages in quartz deposits in joint-like fractures in nearby Cambrian Flathead Sandstone compared with temperatures inferred from burial history are compatible with north- and NW-striking fractures in those rocks forming prior to and during emplacement of the Laramide Buck Mountain reverse fault (Forstner et al. 2018), so similarly orientated fractures in the underlying granite may have formed contemporaneously.

Regular spacing typifies many fracture arrays (e.g. Narr & Suppe 1991; Bai et al. 2000) and regular spacing is sometimes a useful assumption for practical applications (e.g. Wehunt et al. 2017). But irregular to highly clustered fracture patterns are increasingly recognized as important for subsurface fractures (Questiaux et al. 2010; Li et al. 2018). Clustered fractures are viewed as a key to the hydromechanical properties of crystalline rocks (Gabrielsen & Braathen 2014; Torabi et al. 2018). Clustered patterns typify joint patterns in many granites (Genter & Castaing 1997; Ehlen 2000). In granite, both closely spaced joint arrays and patterns arising from the development of wing-cracks are known (Martel 1990; Mazurek et al. 2003). Documenting and quantifying such clustered patterns is a useful objective.

Our exposures may be an analogue for fracture zones in subsurface granites where the effects of fractures on fluid flow and rock strength are issues. For example, the Lancaster reservoir off Shetland produces oil from fractured granite, gneiss and other high-grade metamorphic rocks (Slightam 2014; Belaidi et al. 2018). Basement rock types here closely resemble Wyoming Province crystalline rocks. Moreover, the size and shape of the off-Shetland fractured mass is comparable to that of the Teton Range (cf. Slightam 2014; and Fig. 1). For the off-Shetland example, seismic data supplemented with outcrop observations reveal fracture zones in a range of orientations. Some are 15–75 m wide, comparable to the cluster widths we found. Our example, or the NCC method applied to NW Scotland analogues, could provide insight into the hierarchical arrangement of fracture clusters that may exist at and below the resolution limits of seismic methods.

Microfracture populations are present in many granites (Anders et al. 2014). In some granites, a wide range of lengths can be described with power laws (Bertrand et al. 2015). We cannot rule out that microfractures are present in our outcrop, although the fracture set we measured appears to have a narrow aperture size range. Indeed, the narrow aperture size of visible fractures in this set falls within some definitions of microfractures, although the generally long trace lengths do not. In our outcrop, the lack of staining away from large joints suggests that little fluid flow has occurred through microfractures (those small enough to require microscopy to detect: Anders et al. 2014), a pattern that contrasts with the Lancaster Field, where microfractures are conduits for flow and storage (Belaidi et al. 2018).

For assessing fracture networks, 1D observations are frequently all that can be obtained from the subsurface, and 1D data are useful elements in pattern assessment, and as a link between subsurface observations and the 2D and 3D descriptions and topological analysis of geometrical, as well as spatial, characteristics of networks, including connectivity (Sanderson & Nixon 2015) or to seismic observations (Hu et al. 2018). Our scanline site is notable for providing a clear and readily accessible view of joint clustering in granitic rocks, and as a demonstration of the NCC method. Fracture patterns in this exposure and in other granites have a considerable 2D and 3D complexity over a range of scales (Le Garzic et al. 2011; Bertrand et al. 2015), and so our 1D scanline results are only a first step towards 3D descriptions where cluster dimensions and shapes can be rigorously defined (e.g. discussion in Marrett et al. 2018). Statistically defined cluster boundaries, such as those noted in Figure 8, provide non-arbitrary criteria for delineating dimensions of clusters (swarms, corridors).

Self-organized clustering over a length scale range from 10⁻² to 10¹ m is indicated for north-striking opening-mode fractures (joints) in Late Archean Mount Owen Quartz Monzonite by analysis with the recently developed CorrCount software and the NCC method. These are regional joint patterns in a structural setting >2 km from folds and faults. Using a protocol for interpreting NCC and intensity diagrams, we find that joint clusters are a few to tens of metres wide, and are separated by areas as much as 70 m wide that are largely devoid of fractures. Within large clusters are smaller self-similar spacing patterns.

By normalizing frequencies against the expected frequency for a randomized sequence of the same spacings at each length scale, the degree of clustering and the margins of clusters (or corridors) of anomalously closely spaced, narrow arrays are rigorously defined. Compared with a few clustered patterns in sedimentary rocks that have been described with the NCC method, out example is markedly more clustered, as shown by high exponents and steeper slopes above the 95% confidence limits. Our pattern may reflect the propagation of corridor-like process-zone joint clusters by the mechanism described by Olson (2004) and local wing-crack patterns forced near the tips of pre-existing fractures, a process common in other joint arrays in granite (e.g. Martel 1990).

Our NCC analysis of spacing datasets from other fractured granites, based on measurements previously published by Ehlen (2000), finds patterns that range from indistinguishable from random to marked hierarchical clustering, similar to those in the Wyoming example we report. These examples illustrate the utility of the NCC method over other approaches for describing the spatial arrangement, quantifying clustering and identifying patterns that are statistically more clustered than random.

We are grateful for research access to the Jedidiah Smith Wilderness to Diane Wheeler and Jay Pence, the United States Forest Service and to Grand Teton National Park for permission under National Park Service permits GRTE-2014-SCI-0021 and GRTE-2017-SCI-0072. Stephanie R. Forstner and Finn Tierney assisted with fieldwork. I.A.M. Laubach, E.C. Laubach and A.M. Laubach provided essential field support. We appreciate anonymous reviewer comments on the manuscript, comments on the research from L. Gomez and J.N. Hooker, and insights from members of the 2018 trans-Tetons field trip, including Sergio Sarmiento, Tiago Miranda, Thiago Falcão, Marcus Ebner, Rodrigo Correia and Lane Boyer.

Our research on fracture pattern development is funded by grant DE-FG02-03ER15430 from Chemical Sciences, Geosciences and Biosciences Division, Office of Basic Energy Sciences, Office of Science, United States Department of Energy, and by the Fracture Research and Application Consortium.

Correction notice: the ORCID for Qiqi Wang was added.