En echelon fractures and veins are among the most common and distinctive geological structures, yet their three-dimensional forms and relationships to surrounding structures are commonly unclear. X-ray computed tomography (CT) offers an unrivaled ability to examine structures within rocks in three dimensions, and it is applied here to a sample of drill core from the Marcellus Shale of southwestern Pennsylvania (USA). CT images yield qualitative and quantitative data on the transition from a pyrite-rich planar vein to an en echelon veinlet array, and on the heterogeneity of veinlets within the array. Using a combination of three- and two-dimensional images, geometric data, and traditional petrography, we identify a range of veinlet shapes consistent with deformation during formation of an antitaxial graphite-calcite-pyrite vein system. Each of the veinlets is rooted in the underlying planar vein where it is narrowest. The transition from planar vein to en echelon array coincides with a change in bedding, suggesting that competency contrasts between adjacent beds controlled the fracture morphology. Veinlets initiated as short, lenticular fractures at ∼45° to the planar vein before lengthening, dilating, and rotating. None of the veinlets are strongly sigmoidal, nor is there measurable offset across the margins of the planar vein; therefore, finite non-shear strain was very limited, and fluid overpressure–induced fracturing during burial and diagenesis is probably the most likely process for fracturing and vein formation.


Arrays of en echelon fractures and veins (“tension gashes”) are ubiquitous structures in many geological settings and are exhibited at millimeter to kilometer scales (Pollard et al., 1982). The origins of en echelon fractures, which are among the most easily recognized structures in geology, are debated, in part because of a lack of understanding of their three-dimensional (3-D) geometry. The potential for dilation across en echelon fractures makes them important for understanding syn-deformation fluid flow, pressure solution, mineralization, and vein formation. They are usually observed on single surfaces (e.g., Beach, 1975; Sanderson and Peacock, 2019), therefore understanding of their geometries into and out of the plane is commonly limited to outcrops with multiple oblique surfaces (Smith, 1995), hand samples and loose blocks (e.g., Nicholson and Ejiofor, 1987; Thiele et al., 2015), and analog models (e.g., Mandal, 1995). This paper documents an enigmatic, small en echelon array of millimeter-sized, weakly sigmoidal veins in a sample of shale using 3-D X-ray computed tomography (CT), where, for only the second time (cf. Craddock and van der Pluijm, 1988), we capture the transition from a single planar vein to an en echelon vein array (cf. Lutton, 1971).

En Echelon Arrays

En echelon fractures and veins form during tectonic or locally induced brittle-ductile deformation (Fossen, 2016). Sustained deformation may generate multiple generations of the superimposed en echelon fractures (e.g., Ramsay and Huber, 1987) or destroy the host rock’s coherence to form a fault breccia (Mandal, 1995); conjugate pairs of arrays are very common (e.g., Roering, 1968; Beach, 1975; Ramsay and Huber, 1987). The transition from simple, planar fractures to arrays of en echelon fractures is obscured by the two-dimensional (2-D) nature of most exposures. Lutton (1971) drew an analogy between hackle fringes associated with joint propagation surfaces and en echelon fractures, thereby implying that en echelon arrays were the distal margins of an area of arrested joints (Pollard et al., 1982). If this is correct, the transition from planar fracture to en echelon array may exhibit other features characteristic of tensile joints: plumose structures, hackle lines, and ribs (e.g., Fossen, 2016, his figure 8.17).

Arrays of subsidiary fractures and veins oblique to a planar discontinuity (i.e., joint, fault, ductile shear zone) may form in very different stress regimes that are difficult to distinguish. The shapes of individual fractures or veins inform on the orientation of the strain axes during deformation; however, whether or not their occurrence in en echelon arrays is significant is still debated (see Olson and Pollard, 1991). Subsidiary fractures and veins formed in coaxial deformation form oblique-acutely (<45°) to the planar discontinuity and straight (“wing cracks”; e.g., Willemse and Pollard, 1998; Mutlu and Pollard, 2008), and become longer and more curved as deformation continues (Lee et al., 2016). Wing cracks usually only occur on one side of a discontinuity and are commonly irregularly spaced along its length.

Subsidiary fractures and veins formed in non-coaxial deformation initiate oblique (∼45°) to the shear zone boundaries (Lajtai, 1969; Durney and Ramsey, 1973), and are then rotated toward parallelism (Shainin, 1950; Fossen, 2016) and deformed to form sigmoidal fractures or veins. Alternatively, en echelon arrays may form by the nucleation of grain-scale heterogeneities into microcracks and eventually into macroscopic fractures (Olson and Pollard, 1991), weakening the rock and subsequently localizing shear strain. Many centimeter- to meter-scale en echelon veins are characteristically sigmoidal rather than planar (Ramsay and Huber, 1987), suggesting that brittle fracturing is followed by progressive folding of initially planar fractures within a broader, non-coaxial shear zone (e.g., Cloos, 1955; Passchier and Trouw, 1996). Such a scenario might occur at low but sustained strain rates allowing ductile deformation following instantaneous brittle failure (e.g., Beach, 1980; Rickard and Rixon, 1983), or where a shear zone crosses lithological boundaries between rocks with different strengths (i.e., competency contrasts).

The fractures are thought to initially open (i.e., dilate) in the σ3 direction and lengthen by tip propagation in the σ1 direction (e.g., Beach, 1975; Nicholson and Pollard, 1985; Passchier and Trouw, 1996; Fossen, 2016). How fractures lengthen in the σ2 direction is uncertain. Similarly, how fractures may or may not have connectivity with each other is uncertain (e.g., Thiele et al., 2015).

During progressive deformation, the central parts of the fractures (e.g., Beach, 1975) and the intervening rock masses (Nicholson and Pollard, 1985; Becker and Gross, 1999) are rotated from the σ1 direction to become oblique (Passchier and Trouw, 1996); however, the fracture tips remain parallel to σ1 (Fossen, 2016, his figure 8.30), hence the fractures become sigmoidal. In the “passive deformation” model (Ramsay and Graham, 1970; Hancock, 1972; Durney and Ramsey, 1973; Beach, 1975), the fractures deform in response to the strain gradient across the shear zone. In the “rock bridge” model, the fractures and intervening rock masses shorten and buckle as they rotate (Nicholson and Pollard, 1985; Nicholson and Ejiofor, 1987; Nicholson, 2000), potentially changing the magnitude of dilation across the fractures; eventually the “rock bridges” cannot shorten further and the shear zone hardens, perhaps initiating a conjugate array or the development of a new set of tension gashes (see discussion of the different models in Lisle [2013]).

Alternatively, some numerical and analog models suggest that sigmoidal fractures may initiate with that shape due to inhomogeneous strain between closely spaced fractures during coaxial or non-coaxial deformation (“mechanical fracture interaction”; Olson and Pollard, 1991; Mandal, 1995). A competition may develop between the local stresses around crack tips and the external stress field (Thomas and Pollard, 1993). Therefore, ductile deformation may be of little significance in many arrays of sigmoidal veins, especially weakly sigmoidal ones.

Therefore, understanding of how the shapes, sizes, and spacings between fractures evolve throughout deformation as well as throughout a 3-D shear zone is incomplete and is necessary for future investigations. Herein, we apply 3-D X-ray scanning of an array of pyrite-rich, en echelon veinlets and the planar vein that they join with to investigate the 3-D morphology and vein architecture.

Case Study: Coldstream 1MH Well

The Marcellus Shale gas play in the upper Paleozoic Appalachian Basin of the eastern United States (e.g., Zagorski et al., 2017) has experienced rapid development coinciding with a revolution in drilling and rock-characterization technologies. As a result, the Marcellus Shale can now be drilled and exploited rapidly, and whole-core sections can be extracted and analyzed within days of drilling. Successful exploitation of shale gas plays requires a detailed, 3-D understanding of the geometries, scales, and networking of preexisting fracture sets (e.g., Gale et al., 2007, 2014). Therefore, whole-core CT is routinely applied to describe (e.g., Wildenschild and Sheppard, 2013) and quantify fracture sets in rocks and in the products of experimentally induced fracturing episodes.

Our sample of Marcellus Shale is from the vertical Coldstream 1MH well (American Petroleum Institute well number 37-033-26848) in Clearfield County, central Pennsylvania (41.112667°N, 78.406590°W), where the Marcellus Shale is at 2134.21–2187.24 m below surface. The well was drilled in 2011 and remained an active gas producer as of November 2019. The core was analyzed at the National Energy Technology Laboratory (NETL) in Morgantown, West Virginia, and drilled by Energy Corporation of America (now Greylock Energy).

A characterization of the full Coldstream 1MH core (Crandall et al., 2018), including whole-rock X-ray fluorescence data, petrophysical data, and a complete archive of CT scans, is available at https://edx.netl.doe.gov/group/core-characterization. The industrial-CT (indCT) scans of specific core sections used in this study are available in the same digital archive. The core is cataloged in the System for Earth Science Sample Registration (SESAR; IEDA, 2018), and each 1 m section is assigned a unique International Geo Sample Number (IGSN).

The Middle Devonian Marcellus Shale is an organic-rich shale interbedded with fossiliferous calcareous shale and silty mudstone (Lash and Engelder, 2011). Fine to medium silt-sized quartz and mica grains are present throughout, with patchy calcite cement and disseminated pyrite and organic material (e.g., Lash and Blood, 2014). The shale is pervasively fractured (e.g., Wilkins et al., 2014), with calcite- and pyrite-filled veins both parallel to and discordant with bedding. The presence of pyrite along vein walls enhances the ability of CT to accurately display their morphology. This study focuses on one such bedding-discordant vein.


CT Methodology and Image Processing

The Coldstream 1MH core was scanned by a medical-CT scanner at NETL Morgantown (Toshiba Aquilion TSX-101A/R; Crandall et al., 2018) with a pixel resolution of 430 × 430 × 500 μm. Higher-resolution scans were performed on specific subsections of the core with the NorthStar Imaging Inc. (NSI) M-5000 Industrial-CT system at resolutions of 29.7–70.6 µm3. The indCT scans were performed at a voltage of 185 kV and a current of 400 µA, with scans consisting of 1440 radiographs that were reconstructed into volumes using NSI’s EFX reconstruction software and exported as TIFF stacks. The indCT scans required the core to be scanned in two separate sections, which were then combined using ImageJ software (Rueden et al., 2017). Radiographs were processed in Northstar EFX-CT software to generate grayscale TIFF stacks ready for image processing and analysis.

TIFF stacks were analyzed using the freely available software ImageJ2/Fiji (Schmid et al., 2010; Schindelin et al., 2012). In CT images, the grayscale is primarily controlled by the atomic mass of the material. Reconstructed digital models can be viewed as a series of 2-D slices in any cardinal orientation or as fully rotatable, digital 3-D models. Montages of 2-D slices depict the fracture within the surrounding rock mass, revealing the relationship between bedding, veins (annealed fractures), and open fractures. Three-dimensional models are presented as animations of the segmented vein where the surrounding rock mass has been removed (segmented) in ImageJ2/Fiji.

Structure-from-Motion Photogrammetry

The core sample was photographed using a Sony Cyber-Shot DSC-HX80 18.2 megapixel digital camera mounted to a tripod. The sample was staged on a turntable within a light box illuminated from above, front right, and front left. Fifty-eight (58) and 63 images were taken in three orbits of the top and bottom of the sample, respectively. The images for each half were reconstructed as point clouds and meshes using structure-from-motion photogrammetry (Favalli et al., 2012; Westoby et al., 2012) in Agisoft Metashape software, before being cleaned of artifacts and merged to form a single 3-D model.



Our sample (Fig. 1) is an ∼50-mm-thick piece of a two-thirds slabbed 101-mm-diameter core from ∼2166 m below surface in Coldstream 1MH core 2, box 29 of 31 (IGSN IENTL00VU; https://app.geosamples.org/sample/igsn/IENTL00VU). The upper surface exposing the en echelon vein array (Fig. 1A) is bedding-parallel but also the top of a new core section (Crandall et al., 2018, their figure 26); therefore, continuity with the preceding core is uncertain. The lower surface (Fig. 1B) is also bedding-parallel and appears to have formed along a nonmineralized (i.e., open) fracture probably formed during coring or core extraction. Similar open fractures are pervasive throughout the sample (Fig. 2). The adjacent core pieces do not host any significant veins (Crandall et al., 2018, figure 26), so again, continuity within the core is uncertain.

The sample is composed of carbonaceous shale with interstitial calcite, pyrite nodules, and small veins. Diffuse bedding is picked out in the indCT images by subtle changes in brightness (Fig. 2; Video 1) probably related to different calcite contents. Bedding is laterally continuous across the width of the core, and the beds are not visibly offset across the width of a bisecting vein (i.e., the vein is a filled fracture, not a fault).

Vein Morphology and Petrography

The vein transects the width of the core sample and therefore has a minimum length of ∼10 cm without its ends being preserved. The vein is divided into two distinct parts: (1) a single, planar vein exposed at the base of the sample (Figs. 1B, 3A–3B), and (2) a complex, en echelon vein array at the top (Figs. 1A, 3C).

Planar Vein

The single, planar vein is observed to narrow and bifurcate upwards into the en echelon array, where it is exposed on the outer edges of the upper 2 cm of the sample (Fig. 1C; Video 1). The single, planar vein also exhibits a prominent, laterally continuous bend in the upper half of the sample (Figs. 1C, 2). Video 1 shows that on one side of the sample, the bend is adjacent to a pyrite nodule, and that on the other side, it nucleates a minor, bifurcated en echelon vein segment that reaches the upper surface. Figure 2 shows that in the center of the sample, the bend nucleates a thinner, blind, splaying vein that reconnects with the main vein.

The planar vein is antitaxial and composed of three minerals distributed in three semi-parallel zones: margins—pyrite and calcite; mid-vein—calcite with pyrite; and central lens—graphite (Figs. 3A–3B). The pyrite occurs as discrete, equant, ∼200–500-µm-diameter grains best observed in the segmented indCT data (Videos 2 and 3). The blocky pyrite forms a network of partly interconnected grains on each side of the vein (Fig. 3B; Video 2) and also occurs as disseminated grains within calcite in the vein’s interior. Elongated blocky calcite occurs throughout the vein (Fig. 3B): interspersed with pyrite at the margins, and as the overwhelmingly dominant phase in the main part of the vein in the lower half of the sample where the vein is widest. As the vein narrows upwards, pyrite appears to become equal volumetrically to calcite (e.g., Fig. 2). At the transition from planar vein to vein array, pyrite is very sparse along the center of the vein but abundant at either end (Video 2). Graphite is interspersed with calcite and occurs in an ∼40-mm-long, ∼0.5-mm-wide central lens within the vein (Fig. 3A). Like the calcite, graphite is not discernible in the indCT data due to its low density. A detailed discussion of the paragenesis of the vein is beyond the scope of this study; however, the lenticular shape of the graphite suggests that it is the oldest part of the antitaxial vein (e.g., Passchier and Trouw, 1996, their figure 6.12).

There is a weak preferred mineral alignment within the calcite, and possibly the pyrite, oblique to the margins of the planar vein (Fig. 3D). The alignment is picked out by the prominent cleavage faces of calcite twins (e.g., Burkhard, 1993) ∼50° from the margins of the vein. Some pyrite crystals adjacent to calcite have faces ∼50° from the margins of the vein, producing an asymmetrical sawtooth pattern.

En Echelon Vein Array

The en echelon vein array (Fig. 3C) is composed of ∼20 veinlets of different shapes and sizes, not all of which are exposed on the sample’s upper surface but which are all visible in the indCT data (e.g., Fig. 4; Video 2). The array is ∼1 cm wide and is parallel to the single, planar vein below it. Sixteen discrete en echelon veinlets were identified on the upper surface (Fig. 3C) and in multiple image slices within the 2-D TIFF stack, and measured. The shapes and sizes of en echelon veinlets are described below.

Most of the veinlets in the en echelon array are single segments with varying degrees of curvature. Veinlet 12 is composed of at least three interconnected smaller segments that produce a complex sigmoidal shape (Fig. 3E). The tips of each veinlet are typically gently curved and asymptotic with the array margin and come to a clearly defined point (Fig. 3F). Several of the veinlets are composed of linked vein segments and intravein fragments (e.g., Becker and Gross, 1999), and have characteristically jagged shapes, for example vein 12 (Fig. 3E).

The rock masses between veinlets have sharply defined edges, show no internal disruption, and vary in thickness (Fig. 3G). They are not disrupted by fractures or pressure-solution seams linking veinlets (cf. Peacock and Sanderson, 1995; Thiele et al., 2015) nor are they obviously folded, but there is an absence of visible or CT-imaged passive markers in an appropriate orientation.

The en echelon veinlets are composed of pyrite and subsidiary calcite (Fig. 3F); no graphite has been observed. The pyrite is continuous along and across each veinlet (Videos 2 and 3), rather than as discrete grains on either side as in the planar vein (cf. Fig. 3D), with the widest grains being 100–200 µm across. Disseminated pyrite grains are notably absent.

Quantitative Data

Measurements of the veins within the array were completed using ImageJ2/Fiji’s built-in measurement toolset. The 12 numbered veins (Fig. 3C) were measured in regularly spaced slices to track changes with depth into the sample. Measurements include (Fig. 5): vein or vein-array width, the point-to-point length between veinlet terminations (A), the axial length between veinlet terminations (B), the distance between adjacent veinlets along the projected medial axis of the array (Z), and the angle counterclockwise between the axial trace of the veinlet and the medial axis of the array (φ). α is the angle between the projections of the axial trace at the midpoint and the termination of the veinlet (Fig. 5); it is 180° in lenticular veinlets and <180° in sigmoidal veinlets.

Veinlet orientation, geometry, and scale data are presented in a series of cross-plots (Figs 610) as described below.

Vein and Vein-Array Width

The single, planar vein at the base of the sample (Fig. 3A) is 2–2.3 mm wide (Fig. 6). The vein narrows steadily through the sample to a minimum of 0.4–0.7 mm, at a rate of 0.06–0.1 mm/mm. At −6.5 and −5 mm, respectively, the two measured parts of the vein widen steadily and rapidly (1.4–1.7 mm/mm) into the vein array. The maximum width of the en echelon array measured is ∼9.5 mm.

The pattern is remarkably similar in both measured sections (Fig. 6), and agrees very well with measurements of the vein and vein-array widths on the outside of the sample. The thinnest parts of the vein (−6.5 mm and −5 mm depth) are easily observed in Figure 4 (slices 157–196) and Video 2 where pyrite grains become sparse and the vein is difficult to resolve.

Veinlet Orientation

The ∼10-mm-wide en echelon array is composed of ∼20 discrete veinlets as much as 6.5 mm deep (depth relative to upper surface) (Fig. 3C; Video 2). Orientation data are presented for all 16 measured. The angle φ describes the orientation of each veinlet about the medial axis of the array (Fig. 5). Where φ = 90°, the veinlet is perpendicular to the array. Angle φ ranges from 16° to 58°, with a mean of 36° ± 6° (1 standard deviation) from 124 measurements on 16 different veinlets (Fig. 7). Each veinlet shows a range of φ values with depth, usually a range of ∼10° but as much as 18° (veinlet 1) or as little as 5° (veinlets 9 and 12). Angle φ typically does not vary in response to change in veinlet length (A; Fig. 7A) or veinlet curvature (B/A; Fig. 7B); instead, each veinlet forms a discrete cluster of data. Figure 8 shows geometric data specific to veinlets 1, 7, 11, 13, and 17. Angle φ shows a weak, systematic decrease in the deepest and shortest parts of veinlet 1 (Fig. 8B), and variation in veinlet 13 where φ varies from 45° to ≤35° and back to 45° over its 10 mm extent (Fig. 8K).

Veinlet Length

Veinlets vary in length A (Fig. 5) considerably, from 1200 to 10,500 µm, with most individual veinlets varying (with depth from the sample upper surface) by <2000 µm (Fig. 7A). Mean length A is 5435 ± 2116 µm. Generally, veinlet length decreases systematically with increasing depth; however, veinlet 7 widens before then narrowing (i.e., is widest in the middle; Fig. 8D). Veinlets 1, 7, and 13 vary gradually over a range of 1000–2000 µm; in contrast, veinlets 11 and 17 both have much shorter lengths at their bases than in the body of the veinlet (Figs. 8G, 8M).

Veinlet Shape

The degree to which veinlets are sigmoidal is analyzed in two ways: (1) by comparing the point-to-point length (A) and the axial length (B), and (2) by examining the angle α (Fig. 5). The B/A ratio is a useful shape factor to assess how non-planar (i.e., how sigmoidal) a veinlet is, where the more sigmoidal the veinlet is, the higher the value. All measured veinlets have B/A values <1.1 (Fig. 7B) except for a single measurement on veinlet 4 (B/A = 1.16), and most are <1.05. Mean B/A is 1.02 ± 0.02; therefore, the veinlets are weakly sigmoidal. Angle α varies from 140° to 180°; with a mean of 160° ± 8° (Fig. 9A); therefore, again the veinlets are weakly sigmoidal. Generally, α decreases from 180° and B/A increases, i.e., curvature increases, as φ increases (Fig. 7B); however, given the limited range of both B/A and α values, a statistically significant correlation is not discernable.

Distance between Adjacent Veinlets

The distance between veins (Z; Fig. 5) is measured between the midpoints along the projected medial axis of the array. Z varies considerably from 25 to 6200 µm, with a mean of 1990 ± 131 µm (Fig. 10). However, with the exception of veinlet 9, the longest and most variable in length, all veinlets have small ranges of Z (range of 500–1000 µm). Shorter veins typically have shorter distances between them (Fig. 10), but otherwise, Z does not correlate with other measured variables.

2-D and 3-D Visualizations

Computed tomography data provide an opportunity to examine geological structures in an engaging and nondestructive way. We take advantage of this potential to generate partial cutaways (Video 1; Video S12), animated segmentations (Videos 2 and 3), and 2-D fly-throughs (Video 4; Video S23). These videos show the 3-D shape, size, and orientation of the veinlets in the en echelon array, which would otherwise not be visible. Video 4 and Video S2 are high- and medium-resolution, respectively, fly-throughs of the sample, parallel to the vein and vein array. Both show the narrowing of the planar vein, the prominent bend in the planar vein, and the bifurcation of that vein to form the en echelon veinlets. Video 4 and Video S2 confirm that the planar vein and the en echelon vein array are contiguous.


Transition from Planar Vein to En Echelon Vein Array

The sample images and indCT data show that the planar vein and en echelon array are connected throughout the sample (e.g., Fig. 1; Video 2). The inflection points in Figure 6 show the depth of the transition from the almost pinched-out single vein to the widening en echelon vein array at −5 mm to −6.5 mm. The pinched zone and the transition to the bifurcated veinlets are evident in Figures 1C and 2, and Videos 1–3. This, coupled with the change in X-ray attenuation in the adjacent rock masses, suggests that a bedding plane within the shale juxtaposes two layers of different competencies (calcite rich versus calcite poor), which responded differently to the propagation of the same fracture.

The cross-sectional morphology of the vein is similar to that of tensile fractures translating toward hackles at the fracture tip (e.g., Lutton, 1971; Pollard et al., 1982). In this scenario, the prominent bend in the vein beneath where it nearly pinches out is a “rib” or “arrest line” on the fracture plane. However, we find no evidence of hackle lines or plumose structures on the planar vein surfaces, although whether they would be coated by pyrite is uncertain. If this analogy is correct, then the fracturing is purely tensile (mode I) and likely to be the result of coaxial shear strain.

Veinlet Heterogeneity and Possible Fracturing Mechanisms

Veinlets in the en echelon array are different from one another in length (A; e.g., Fig. 7A) and distance apart (Z; e.g., Fig. 10), and to a lesser degree, geometry. There is almost one order of magnitude difference in length A between the smallest and largest veinlets. Moreover, while most veinlets steadily decrease in length downwards toward their intersection with the planar vein (e.g., veinlet 13; Fig. 8), some (e.g., veinlet 7) first increase and then decrease, and others decrease unsteadily (e.g., veinlet 17). Although plotting in narrow ranges throughout the entire data set, angles φ and α and ratio B/A, vary nonsystematically and subtly within each veinlet.

The heterogeneity within and between veinlets in terms of size and shape suggests that, at least at the small scales analyzed, veinlets have complex morphologies that change in three dimensions, and that they should not be considered perfectly uniform: neither simple planes nor sigmoids. This agrees with the few 3-D studies of macroscopic en echelon arrays available (Nicholson and Ejiofor, 1987; Thiele et al., 2015), although our examples are less sigmoidal.

The weakly sigmoidal shapes of the veinlets visually match those produced in numerical simulations by Lisle (2013) for non-coaxial shear strain (γ) ≤2.5 and a heterogenous shear zone with a shear profile with high gradients at the margins. The angle φ averages 36° (Fig. 7A), suggesting ∼9° of rotation for most veinlets from their starting orientation of φ ≅ 45° (Durney and Ramsey 1973; Fossen, 2016), assuming that they developed parallel to σ1. This equates to γ = 0.16. However, veinlets 10, 13, and 17 show less (<5°) rotation, and veinlet 1 shows significantly more (∼20°). This could be explained as different magnitudes of finite stress where veinlet 1 is the oldest and has accumulated the most shear strain, and veinlets 10, 13, and 17 are the youngest and least strained. This sequential development could be tested further using isotopic geochemistry and fluid-inclusion thermometry.

However, the lack of visible offset within the shale (Video 4) and the weak sigmoidal shapes of the veinlets (Video 3) suggest that non-coaxial shear strain (i.e., mode II) may not be the primary formation mechanism. Angle α and ratio B/A are both very consistent and low (Figs. 7B and 9), reflecting that none of the veinlets are strongly curved; the only outlier is veinlet 4, where α = ∼180° (i.e., straight) and B/A varies from 1.01 to 1.16, the largest range of B/A values. However, veinlet 4 is the shortest (<2000 µm; Fig. 7A), so small differences in A and B may be exacerbated. The veinlets are similar in shape to many wing cracks formed by fluid expulsion (e.g., Mutlu and Pollard, 2008; Kobchenko et al., 2014; Ma et al., 2016) and modeled by Willemse and Pollard (1998), Lee et al. (2016), and others. This suggests that fluid overpressure during compaction (vertical uniaxial stress) and diagenesis may have been sufficient to fracture already cemented shale, and that the fluids were probably internally sourced. In this scenario, the upward-propagating fracture arrested against a competency contrast boundary within the shale and spawned an array multiple, similar wing cracks (e.g., Craddock and van der Pluijm, 1988). If this is correct, it represents an unusual and, to the best of our knowledge, unique example of a well-organized array of wing cracks propagating away from one another to produce an array difficult to distinguish from a “classic” set of sigmoidal tension gashes.

In conclusion, it is impossible to completely differentiate coaxial versus non-coaxial shear strain from the single small sample examined. The difficulty in distinguishing between coaxial and non-coaxial shear strain in forming en echelon veins is discussed by Olson and Pollard (1991), who noted that both field and laboratory investigations produced similar structures by different mechanisms. This topic remains fertile ground for further investigation, but based on the weight of evidence available, we tentatively favor formation by coaxial shear strain.

Interaction between Veinlets

The spacing between veinlets is positively correlated with veinlet length (Fig. 10). Spacing does not appear to be correlated to the degree of curvature or the angle φ. Spacing therefore seems to be intrinsic to how the fractures initiated, rather than how they evolved. We see no evidence for the large, broken rock bridges between veinlets that would be expected at high magnitudes of finite strain (e.g., Nicholson and Pollard, 1985). Several other studies of en echelon vein arrays have highlighted the importance of pressure-solution surfaces and minor fractures that cross rock bridges (Thiele, et al., 2015) and carbonates (Peacock and Sanderson, 1995), but we see no similar structures (e.g., stylolites) in our data. We suggest that this is because our veinlet array is so small and still connected to the “parent” vein, such that separate discrete fluid pathways did not need to develop between veinlets.


We have analyzed the morphology and structure of a simple planar vein transitioning into a complex, en echelon veinlet array within a single piece of core using nondestructive CT imaging. The parent graphite-calcite-pyrite vein is antitaxial, and it narrows to almost pinching out before bifurcating to form the en echelon array at a bedding plane, possibly due to a competency contrast between beds. The veinlets occur over a range of sizes and spacings, although most are of similar shape and orientation, suggesting weak strain possibly induced by excess hydrostatic pressure, although a far-field tectonic stress component cannot be ruled out entirely. Moreover, the absence of fractures or pressure-solution surfaces between veins indicates that fluid transmissivity was via the adjoining planar vein structure, and that the intervening rock masses deformed ductilely.


This work was completed at the National Energy Technology Laboratory (NETL) with support from U.S. Department of Energy’s (DOE) Office of Fossil Energy Oil & Gas Program. Photographs in Figures 1A–1C are by Karl Jarvis (NETL). CT scans were performed by Bryan Tennant and Scott Workman. This research was supported in part by appointments from the NETL Research Participation Program, sponsored by the U.S. DOE and administered by the Oak Ridge Institute for Science and Education and the West Virginia University Libraries Open Access Author Fund. Reviews by Julia Gale, Nick Hayman, and Associate Editor Terry Pavlis greatly improved the manuscript and encouraged us to reassess and broaden our considerations of different stress-strain regimes and fracturing mechanisms.

1Supplemental 3-D files. From Figure 1D and Video 2. Please visit https://doi.org/10.1130/GES02191.S1 or access the full-text article on www.gsapubs.org to view Supplemental 3-D files. For Figure 1, an annotated, steerable, high-resolution version is also available at Sketchfab, https://skfb.ly/6LYUP.
2Supplemental Video S1. Rotating orthoslice image (64.3 µm/pixel resolution) semi-parallel (long and intermediate axes) to the planar vein. Please visit https://doi.org/10.1130/GES02191.S2 or access the full-text article on www.gsapubs.org to view Supplemental Video S1.
3Supplemental Video S2. 2-D cross-sectional view passing through the sample parallel to the strike of the planar vein (64.3 µm/pixel resolution). Please visit https://doi.org/10.1130/GES02191.S3 or access the full-text article on www.gsapubs.org to view Supplemental Video S2.
Science Editor: Andrea Hampel
Associate Editor: Terry L. Pavlis
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.

Supplementary data