Shale Anisotropy and Natural Hydraulic Fracture Propagation: An Example from the Jurassic (Toarcian) Posidonienschiefer, Germany

Cores recovered from the Jurassic (Toarcian) Posidonienschiefer (Posidonia Shale) in the Lower Saxony Basin, Germany, contain calcite‐filled fractures (veins) at a low angle to bedding. The veins preferentially form where the shale is both organically rich and thermally mature, supporting previous interpretations that the veins formed as hydraulic fractures in response to volumetric expansion of organic material during catagenesis. Despite the presence of hydrocarbons during fracturing, the calcite fill is fibrous, and so, the veins appear to have contained a mineral‐saturated aqueous solution as they formed. The veins also contain myriad host‐rock inclusions having submillimetric spacing. These inclusions are strands of host rock that were entrained as the veins grew by separating the host rock along bedding planes rather than cutting across planes. The veins therefore produce significantly more surface area – by a factor of roughly 5, for the size of veins observed – compared to an inclusion‐free fracture of the same size. Analysis of vein geometry indicates that, with propagation, a fracture surface area increases with fracture length raised to a power between 1 and 2, assuming linear aperture‐length scaling. As such, this type of fracture efficiently dissipates elastic strain energy as it lengthens, stabilizing propagation and precluding dynamic crack growth. The apparent separation of the host rock along bedding planes suggests that the mechanical weakness of bedding planes is the cause of this inherently stable style of propagation.


Introduction
Hydraulic fracturing in shale is important for the sealing capacity of shale layers during natural (Foschi et al., 2018;Petrie et al., 2014) and artificial (Elkhoury et al., 2015) fluid-expulsion events. The mechanical anisotropy of shale is evident from its fissility in hand samples, from abundant examples of layerparallel natural fractures hosted in shales worldwide (Cobbold et al., 2013;Gale et al., 2014;Hooker et al., 2019) and from laboratory experiments indicating that shale has a lower surface energy along bedding planes than across them (Chandler et al., 2016;Luo et al., 2018). This anisotropy could therefore guide fracture propagation into a layer-parallel orientation, improving the effective sealing capacity of shales. This process may be most pronounced for fluid-driven fractures, where high fluid pressures produce low effective stresses and theoretically low differential stresses, due to poroelastic stress coupling (Engelder & Fischer, 1994). Therefore, if imposed differential stresses, for example, from tectonic loads, are small, then the relative importance of rock anisotropy is greater in determining the direction and extent of fracture propagation.
However, in most cases, it is not clear what the past loading conditions were that led to fracturing (Laubach et al., 2019). Subsurface fracturing is theoretically driven by decreases in effective compression (e.g., Gudmundsson, 2011). In shales, such stress changes may be related to tectonic strains (Ghosh et al., 2018;Pireh et al., 2015), fluid overpressures (Meng et al., 2017;Zhang et al., 2016), diagenetic reactions (Hooker, Huggett, et al., 2017), exhumation (Engelder & Gross, 2018), or thermal stresses (English, 2012). Consequently, there is very little constraint on the relative importance of mechanical anisotropy, high fluid pressures, or high horizontal loads for the growth of layer-parallel fractures. Nor is it easy to assess how mechanical anisotropy affects the dynamics of crack propagation, because of poor constraints on how fracture opening and propagation affect the effective stresses in the vicinity of fractures, through their effects on fluid pressure (McDermott et al., 2013;Mohammadnejad & Khoei, 2013;Olson, 2003).
The purpose of this paper is to examine the effects of shale anisotropy on natural hydraulic fractures. We take advantage of the circumstances regarding the formation of a particular set of fractures intersected by cores through the Lower Jurassic (Toarcian) Posidonienschiefer, in the Lower Saxony Basin (LSB), Germany. New observations support a previous study (Jochum et al., 1995) suggesting that the fractures formed in response to fluid overpressures that arose from the maturation of organic matter. It is commonly interpreted that hydraulic fracturing in organic-rich shales forms veins and unmineralized Mode 1 fractures as a direct result of volume change reactions during catagenesis, which have the effect of increasing local fluid pressures (e.g., Lash & Engelder, 2005;Mandl & Harkness, 1987). The four studied cores, from across the Posidonienschiefer in the LSB, are ideally located to test this hypothesis, because they represent a gradation in thermal maturity of the organic matter, from immature to highly mature, with vitrinite reflectance (VR) values from 0.5 to >3. Furthermore, the fibrous mineral fill of the studied veins preserves information about opening kinematics, allowing us to interpret how the veins grew. The present veins dip obliquely with respect to bedding, but fibrous calcite cement deposits show that the veins propagated in en échelon steps, each of which split the rock along bedding planes. Therefore, these natural hydraulic fractures illustrate the stabilizing effect of bedding-plane weakness on hydraulic fracture propagation.

Geologic Setting
The LSB (Figure 1a) is part of the Central European Basin System, which began as a foreland basin to the Variscan orogeny (Betz et al., 1987). The breakup of Pangea in the Late Triassic and Early Jurassic brought rifting and increased accommodation in the LSB (Kockel et al., 1994), primarily via extension along WNWstriking normal faults (Betz et al., 1987). The Posidonienschiefer was deposited, at least in part, during the Toarcian Oceanic Anoxic Event (T-OAE), a period of widespread anoxic-euxinic conditions in marginal  Kockel et al., 1994. (b) Schematic block diagram showing locations of cores in this study and subsurface faults, synthesized from Baldschun et al. (1996). Shaded layer shows approximate maximum paleodepth of the top of the Posidonienschiefer (layer thickness exaggerated). marine environments across both hemispheres (Al-Suwaidi et al., 2016;Dickson et al., 2017;Them et al., 2017) particularly in the northwest continental European Seaway, including the LSB (Jenkyns, 1988(Jenkyns, , 2010. High organic matter sequestration and associated high sedimentary total organic carbon (TOC) values formed in geographically and hydrographically restricted basins, under anoxic-euxinic water column conditions (Jenkyns, 1988).
During the Late Jurassic and Early Cretaceous, relatively rapid subsidence in the LSB resulted in the preservation of a thick, lacustrine-marine succession (Bruns et al., 2013). The Alpine orogeny in the Late Cretaceous reactivated preexisting normal faults in a reverse sense, producing inversion structures, such that exhumation was widespread and to greatest extent within grabens (Bruns et al., 2013).
The localities of the cores used in this study lie within 200 km west of Hanover, Germany ( Figure 1). The Posidonienschiefer in the studied cores lies at present-day depths between 1,200 and 3,500 m (Table 1). Maximum (preinversion) burial depths were estimated from basin modeling by Bruns et al. (2013) and are reported in Table 1. Regional seismic data (Baldschun et al., 1996) imaged numerous faults, including major graben-bounding faults (Figure 1b) that were reactivated during basin inversion.
The Lower Toarcian Posidonienschiefer in the LSB is informally referred to as the Lias Epsilon (e.g., Küspert, 1982), and is underlain and overlain by the Lias Delta (Pliensbachian) and Lias Zeta (Upper Toarcian or Aalenian), respectively (Littke et al., 1991). The Lias Epsilon in the LSB is identified by its elevated TOC and calcite concentrations, such that its exact chronostratigraphic position may vary locally.

Methods
Fractures were identified in the cores as any planar discontinuities containing minerals, slickensides, or organic matter as evidence of formation in the subsurface. Cores tend to break apart along bedding planes, and there are also occasional breaks at a high angle to bedding, often with irregular (nonplanar) attitudes. Breaks of any orientation that lack minerals or slickensides as evidence of formation in the subsurface are attributed to coring or handling and are ignored in this study. Fractures filled or mostly filled by cement are termed veins. Fracture size and position were measured from core photographs, with calibration from physical core observation using a hand lens. For opening-mode fractures, which lack evidence of significant shear displacement, kinematic aperture was recorded in the sense of Marrett et al. (1999) as the total distance between separated fracture walls, regardless of whether the intervening space contains minerals or porosity.
This study focuses on calcite-filled veins, whose fibrous fill (see below) preserves millimeter-scale opening displacements. For such veins, we therefore use the width of calcite cement to reflect the fracture kinematic aperture. This approach may underestimate the true original opening distance of veins if there were any collapse, pressure solution, or other processes that shrunk the veins after their formation. The approach may overestimate the true opening width if expansion during natural exhumation or coring and sampling occurred. However, the core-and petrographic-scale evidence presented below shows mostly intact veins without significant crystal or wall-rock deformation, except that due to the vein opening itself. Therefore, we estimate the error of our size measurements to be ±0.5 mm, based on our image pixel size. Fracture Abbreviation: TOC = total organic carbon. a Based on basin model of Bruns et al. (2013). b Based on seismically detected faults (Baldschun et al., 1996). length was measured where two opposite fracture tips were preserved in the core. Only small shear offsets were noted along fractures, as described below.
Cores were sampled for petrographic observation and geochemical analysis. Petrographic samples were cut into thin sections or polished core slabs. Samples were imaged using an FEI Quanta 650 scanning electron microscope (SEM) at the Department of Earth Sciences, University of Oxford, using a 20-kV beam and 10mm working distance. Crystallographic axes were mapped using an electron backscatter detector attached to the SEM. TOC was quantified using a Rock Eval VI pyrolysis instrument, also at the Department of Earth Sciences, University of Oxford, following the methods in Behar et al. (2001).

Lithology
Cores 1-3 preserve approximately 50 m of gray-brown to black mudrock with intercalated, commonly nodular, limestones ( Figure 2). The bottom 2-5 m of this interval in each core contains unfractured, homogeneous mudrock; above this interval, limestones alternate with shale. The limestones have thicknesses of several centimeters to 1-2 dm, commonly with gradational boundaries and nodular (ellipsoidal) to tabular morphologies. Here, the term shale denotes mudrock having a distinct bedding-parallel fabric of varying clay, calcite, and organic matter composition. The cores tend to break along the sedimentary fabric, and abundant calcite veins lie at low angle to the shale bedding orientation (Figure 3a).
The lowest preserved limestone layers in the studied intervals of each core are coincident with the onset of TOC-enrichment ( Figure 2) at the base of the Posidonienschiefer. Elevated TOC values are generally below 2% in the Lias Delta (below the Posidonienschiefer) and commonly increase to 10-15% in the Lias Epsilon. Within the high-TOC interval, the organic carbon content is generally lower in the limestone intervals than in the shale intervals. Moving upsection, limestones become sparser and are progressively more laminated and less nodular over several meters. Here, laminae denote millimeter-scale layers that are compositionally distinct, primarily by their carbonate mineral content, such that limestone and shale laminae commonly alternate to form centimeter-to decimeter-scale composite sequences ( Figure 3c). Calcareous laminae and TOC both diminish upward, and the change is generally gradual (Figure 2). Core 4 preserves only the section from the top of the Posidonienschiefer down to the laminated interval.

10.1029/2019JB018442
Journal of Geophysical Research: Solid Earth generally framboidal but may be pore-filling or pore-expanding (up to 20%); and apatite (up to 5%). TOC ranges from <1 to 18% ( Figure 2). Laminated and nodular limestone layers contain greater abundances of calcite, both cements and fossils, relative to shale layers, generally at the expense of clay. The top of the organic-rich interval is sharp in Core 1, and more gradational in Cores 2 and 3, but within dark mudrock and not marked by any obvious lithological change in any core. Within the organic-rich section in the latter two cores, carbonate content is reflected in core color, with paler intervals richer in carbonate and poorer in organic matter and clay.
Shale layers within the high-TOC intervals in Cores 2 and 3 contain abundant layer-parallel calcite veins (Table 1). We have no TOC measurements from Core 4, but that core also has abundant veins in equivalent strata, using carbonate beds as correlation markers. Layer-parallel calcite veins are absent in the limestone and calcareous shale layers, and in the low-TOC mudstones of the Lias Delta, underlying the Posidonienschiefer. Shale layers containing layer-parallel veins are not compositionally laminated; however, the long axes of grains in these layers are strongly aligned in a horizontal attitude ( Figure 3).

Calcite Veins
Layer-parallel vein morphologies vary by the angle of the face on which the vein is slabbed (observed). Because the cores are not oriented, the true orientation of any such face is unconstrained. However, individual veins observed on orthogonal cut faces commonly show a characteristic hostrock inclusion pattern on at least one such face. This pattern is a sigmoidal to parallelogram-shaped array of submillimetric vein segments here termed veinlets (Figure 4a-d). The veinlets are thickest, vertically, near the center of the macroscopic vein and thinner toward the tips. Veinlets have close to constant horizontal thickness throughout the vein, and are separated by host-rock inclusions of near-constant, submillimetric thickness, along each inclusion and across the vein.
Vein walls generally make an angle of 10-20°from bedding. Away from the veins, the host-rock grains lie parallel to bedding. Where host-rock is included within the veins, the grains lie parallel to the inclusion (Figure 5a The longest dimension of the fibers is roughly perpendicular to the median line within each veinlet and thus at an oblique angle to bedding. Fibers are rooted at the median lines within each veinlet ( Figure 5). Away from the median lines, crystals become coarser and more equant, particularly toward the narrow ends of the sigmoidal veinlets. Electron backscatter diffraction mapping shows no crystallographic preferred orientation, either among crystals within the median line or among the fibers (Figure 5c).
Vein kinematic aperture ranges from~0.1 to 13 mm (average 2.7 mm, standard deviation 1.8 mm). Veins do not appear to be associated with a significant microscopic-fracture population; no calcite veins were observed to be significantly thinner than the typical thickness of median lines, which is on the order of 0.1 mm (Figure 5b). The maximum tip-to-tip length we observed in core is 72 mm, owing to the 70-mm core width. The maximum aperture observed among small-length veins fully contained within core samples is 4 mm (Figure 6a). Among these small veins, there is a significant, positive correlation between kinematic aperture and length (p < 0.01 based on Pearson's r for linear, exponential, and power-law best fits). The data are insufficient to interpret a linear or sublinear (power law) relationship between aperture and length; neither linear nor sublinear best-fit equations between vein length and aperture can be rejected using chisquared tests.
Consistent with a characteristic veinlet width, the number of veinlets per vein, measured from core observation, varies linearly with total vein length ( Figure 6b). We observe 1.7 veinlets per millimeter of vein length (standard deviation 0.2 veinlets per millimeter). We estimate the uncertainty in veinlet counting at 20%, based on subjectivity of host-rock-strand continuity (Figures 4 and 5).

Journal of Geophysical Research: Solid Earth
Overall layer-parallel vein abundance varies from 0.5 to 19 m −1 within the Lias Epsilon (high TOC) interval (Table 1). Core 1 has the lowest abundance of layer-parallel veins (Table 1). Core 2 has intermediate vein abundance, and Cores 3 and 4 have the highest vein abundance (Table 1). Moreover, analyzing Cores 2 and 3 individually, vein abundance correlates very closely with the TOC in the host rock (Figure 2). Layer-parallel veins are virtually absent away from high-TOC intervals.

Thermal Maturity and Depth of Burial
Rock Eval data ( Figure 2, Table 1) are consistent with previous maximum paleodepth estimates based on VR and other temperature indicators (Bruns et al., 2013) suggesting that thermal maturity varies considerably across the study area. Core 1 is immature, and Cores 2 and 3 are considerably more mature, having hydrogen indices reduced to nearly zero (Table 1). No thermal maturity data are available from Core 4, but its geographic proximity to Core 3 ( Figure 1b) is consistent with a high thermal maturity as well.
The distance between each core site and the closest seismically resolvable fault (Baldschun et al., 1996) is also illustrated in Figure 1b. These were interpreted as inverted normal faults separating regions of strongly disparate maximum burial (Bruns et al., 2013). The degree of exhumation for Cores 2, 3, and 4 is estimated to be on the order of 3,000 m, based on modeled maximum burial depths (Bruns et al., 2013) and the present-day core interval depth (Table 1). These cores were drilled through the same basinscale graben ( Figure 1b); Core 2 was drilled closer to the southern margin, 10 km from the southern fault; Cores 3 and 4 were drilled close to the center, roughly 20 km from either fault. Core 1 is estimated to have undergone a much shallower maximum burial and exhumation of less than 100 m (Table 1); this core was drilled through the footwall of the normal fault bounding this graben to the north, at a distance of roughly 10 km from the fault (Figure 1).

Hydrocarbon Maturation as Vein Formation Mechanism
The observed intensity of veins within the high-TOC interval ( Figure 2, Table 1) is not correlated to fault proximity. In previous studies, fracture intensity has been shown to increase near larger-scales structures such as buttress anticlines (e.g., Glen et al., 2005) and otherwise isolated faults . Such fault-zone fractures are interpreted as a manifestation of progressive brittle strain during fault slip (Ferrill et al., 2019;Laubach et al., 2014). In the present cores, however, the veins are not systematically distributed around the seismically resolvable, grabenbounding faults, nor are there appreciable changes in bedding attitude within the fractured core intervals. We therefore find little evidence that the veins represent fault-zone damage.
Rather, there is a strong correlation between organic matter, thermal maturity, and layer-parallel fracturing (Figure 2, Table 1). Variation in maturity throughout the study area is thought to reflect variation in thermal exposure and is likely linked to structural position: the WNW-striking normal faults produced a horst and graben architecture to the basin prior to inversion during the Late Cretaceous (Bruns et al., 2013), with downthrown fault blocks buried to a significantly greater depth than adjacent blocks. As organic rich rocks mature, the degree of conversion of the original kerogen to

Journal of Geophysical Research: Solid Earth
petroleum-bitumen products also increases. If the expulsion of the petroleum products had been hindered by the generally low permeability of the host source-rock, it is possible that the intensity of the hydraulic fracturing driven by the maturation process would also have increased (see Mandl & Harkness, 1987). A low-permeability, generally closed fluid system is also consistent with previous laboratory studies of shales (Neuzil, 1994), and with the vein fills comprising minerals that are also abundant within the host rock . Our observations therefore support those of Jochum et al. (1995) that layer-parallel veins were formed by volume increase and overpressure related to the maturation of organic matter; i.e., the veins are natural hydrofractures driven by catagenesis.

Growth of Sigmoidal Veinlets 5.2.1. Antitaxial Cementation
In general, veins comprising a median line surrounded by continuous, fibrous fill that extends to the vein walls are interpreted to have filled antitaxially (Bons et al., 2012). That is, such veins originate as thin veins at the median line and then simultaneously widen and fill by adding cement at the vein-wall interface.
Previous numerical models have recreated this antitaxial texture by assuming that crystals at the median line have random crystallographic orientations (Ankit et al., 2015;Hilgers et al., 2001). These models then posit that cementation occurs by epitaxial overgrowth of these randomly oriented seed crystals, such that crystal growth passively fills in the space between the fiber tip and the vein wall. Accordingly, if the walls of a growing vein move apart slowly enough, each crystal will be able to grow fast enough to fill that ephemeral space regardless of crystallographic orientation, and fibrous fill will result. In contrast, if the vein opening rate is faster, only those crystals having fast growth axes at the high angle to the vein wall will keep up with the opening rate, and a coarser fill style will develop, which also has a crystallographically preferred orientation (Ankit et al., 2015;Hilgers et al., 2001;Lander & Laubach, 2015).
The fibrous texture of the vein fill observed here, combined with the lack of crystallographic preferred orientation, therefore supports previous models of cement fill, where cementation was simultaneous with, and at the same rate as, vein opening. The vein opening rate was presumably a function of the vein driving stress, which in this case was the sum of the overburden and the internal hydrocarbon fluid pressure.

Nature of Host-Rock Strands
The veins include myriad strands of host rock. Two competing interpretations of these strands are that they are inclusions, which were incorporated into the veins during vein growth, or, alternatively, that they are the result of fracturing of initially continuous calcite veins. We prefer the former interpretation, for two reasons.
First, the long axes of grains within the host rock, away from fractures, tend to lie parallel to bedding, whereas within host-rock strands, grain long-axes are rotated toward parallel with the strands (Figure 5a, b). Thus, it appears that the sedimentary fabric was rotated within tilted host-rock strands, as opposed to being cleft or truncated by veinlets, or sheared as a result of injection between veinlets. The rotation of the host-rock fabric is consistent with inclusion of host-rock between veinlets that grow against one another while propagating along, and not across, bedding planes.
Second, the host-rock strands delineate the boundaries between calcite crystal domains that are undeformed and have contrasting textures and crystallographic axis orientations. Figure 5c shows that fibers commonly terminate not against the outer vein wall but against host-rock inclusions between the median line and the outer wall. Within the tips of veinlets, outward of host-rock inclusions, the vein fill becomes coarser. If the veins were composed of continuous, fibrous calcite before introduction of host-rock strands, then we would expect general continuity of fibers on either side of those strands. Fibers could be rotated during a hypothetical shearing episode, but we see no significant deviation of the attitude of fibers throughout the veins. Nor The absence of fibers within the narrow ends of the veinlets (Figure 5c) suggests that continued growth (opening and propagation) of the veinlets was filled by cements overgrowing progressively coarser substrates toward the veinlet tips. As with fibers, we might expect to see crystallographically parallel blocky crystals across host-rock strands, if those strands formed by the fracturing of previously continuous calcite crystals. And although some crystals on either side of strands have similar crystallographic orientations (Figure 5c), there are abundant examples of highly contrasting orientations across even relatively thin strands, judging by the inverse pole figure (color key, Figure 5c). This observation favors the interpretation that the host-rock strands were incorporated between originally distinct calcite crystals, rather than having been translated in among previously continuous calcite crystals.
Networks of host-rock strands amid vein-like calcite aggregates are also present in cone-in-cone structure (e.g., Cobbold et al., 2013). Analogous observationsnamely, highly contrasting crystallographic orientations between undeformed crystals separated by thin host-material strandswere used to argue for incorporation of host sediment (Maher et al., 2017) or host rock (Hooker & Cartwright, 2018) during cone-incone growth and not afterward.

Growth Reconstruction
To understand how these veins formed, and the role of shale mechanical anisotropy in their development, we present a growth reconstruction model (Figure 7) that accounts for several key observations. These  Figure 5c). Median line is dashed. Black semicircles show offset of original circles before formation of veinlets. Thin horizontal lines represent bedding planes; regions of closely spaced bedding planes illustrate the original stratigraphic thickness, based on the assumption that thickness is conserved within the (rotated) inclusions in between veinlets. Wide bedding plane spacing is drawn so that both the bedding line-lengths and the total host-rock area are conserved. The true strain distribution around veins is unclear and is likely compensated, at least in part, by neighboring veins. include (i) the antitaxial formation and (ii) syn-growth incorporation of host-rock strands mentioned above, as well as (iii) positively correlated aperture and length (Figure 6), (iv) veinlets that propagate along bedding planes and not across them (Figure 5b), (v) median lines that lie oblique to both bedding and the vein long dimension (Figure 4c,d), and (vi) a characteristic veinlet width, which does not systematically vary with position inside veins or with vein size (Figure 4).
We interpret that veins in the Posidonienschiefer of the LSB initiated as single veinlets, which separated the host rock at a location now marked by the median line of the centermost veinlet (Figure 7). The initial veinlet widened and lengthened along the same bedding plane, filling with fibrous calcite as it grew. But, owing to nonparallelism between the applied stress field and the sedimentary fabric, the vein propagation direction was oblique to bedding. After some enlargement, rather than cutting across bedding planes, a new veinlet (Veinlet 2 in Figure 7) formed near either tip of the original veinlet, beginning an en échelon array.
The vein continued to lengthen in the same manner, with a new veinlet added at either tip once the preceding tip-veinlets reached a certain size. This addition of new veinlets resulted in a regular spacing of veinlets as the vein grew (Figure 7). The earlier-formed veinlets continued to widen as the veins grew; widening was accomplished by progressive separation of the host rock along bedding, i.e., by propagation such that the outer veinlet wall remained parallel to the sedimentary fabric and such that the internal, dipping veinlet wall remained parallel to the sedimentary grain within the rotated host-rock strand. Note that the marker horizons in Figure 7 are continually entrained between the same veinlet pair as the vein grows.
In order to illustrate the veinlet-scale growth of veins in Figure 7, we simplify the median line as a discrete surface through the center of the vein. SEM images (Figure 5a-c) show that the median line is actually a zone of microscopically blocky crystals having myriad tiny host-rock strands included among them. A subset of these inclusions extend to the outer vein wall, forming the veinlet boundaries. Marker horizons in the figure are drawn so as to conserve length, to bend within the vein's host-rock inclusions, and to lie horizontal away from the vein. This approach illustrates the mixed kinematics (opening plus left-lateral shear, as drawn) inferred from the geometry of the crystals and host-rock inclusions. However, the true strain distribution around the veins is unclear, owing in part the to the homogeneity of the host sediment. The true displacement field around the veins is likely compensated, to a large extent, by nearby veins.
We suggest that the veinlets in Figure 7 are longer and blade-shaped in the third dimension, into the page. This would account for the long, continuous shape of vein segments in orthogonal cut faces (Figure 3b, smaller fractures in Figure 3a, 4e-f).
A potential alternative reconstruction model might hold that each macroscopic vein formed by the initial development of an array of uncemented microfractures, which variably grew taller in a second phase of calcite precipitation-related expansion. However, this interpretation is belied by the fact that longer, thicker veins have veinlet widths and spacings similar to those of small veins (Figures 4 and 6b). If this pattern emerged from en échelon arrays of tiny fractures of equal length, with overall array-length a linear function of the number of tiny fractures present, then the final aperture of each calcite vein would seem to have been predestined by the number of tiny fractures in the array. For example, suppose there are two arrays of microfractures alike in all respects except one is formed of 50 fractures and is 30 mm in total length and the other is formed of 100 fractures and is 60 mm in total length. This alternative hypothesis would imply that once these arrays evolve into veins, the latter vein would have twice the final aperture of the former, because it had twice the initial length. In contrast, by our model, each vein grew by widening and adding new veinlets at its tips at regular intervals of growth, accounting for the observation of more veinlets within longer, thicker veins.

Vein Mechanics
The interpretation of vein growth, detailed in the kinematic model outlined above, implies that weakness along the bedding fabric in shales can exert an important control on natural hydraulic fracture propagation. Ideally, opening-mode fractures propagate in a plane perpendicular to the minimum compressive stress direction (e.g., Cosgrove, 1995). Where this direction is misaligned with the pole to the bedding fabric, the fracture propagation path and total distance will reflect an energy minimization between the imposed stress field and the material weakness. Understanding this phenomenon will help to better predict the behavior of induced fractures in the subsurface.

Journal of Geophysical Research: Solid Earth
To quantify the mechanical effects of shale laminae on fracture propagation, we begin by using an approach based on linear elastic fracture mechanics. We posit that the propagation of the macroscopic veins is analogous to that of natural fractures, based on (i) the crack-like geometry of macroscopic veins; (ii) hooking, branching, and en échelon stepping relationships between macroscopic veins; (iii) the presence of median lines; (iv) vein crystal fibers that are mutually parallel among the many veinlets and perpendicular to median lines; and (v) their inferred fluid-overpressure opening mechanism.
Within an elastic rock body, the energy available to propagate a fracture, by creating fracture surface area A, comes from the elastic strain energy U stored within the rock (Gudmundsson, 2011). Neglecting energy losses in the form of heat, sound, or other kinetic forms, fracturing is the transformation of elastic strain energy into surface energy and any associated plastic deformation of the host. Assuming elastic behavior, the change in elastic strain energy as a result of fracture propagation (i.e., an increase in fracture surface area), is given by (1) where G = 2γ and γ is the surface energy of the rock. The factor of two accounts for fracture propagation creating two surfaces. Griffith (1920) derived the magnitude of U 0 , the amount of elastic strain energy diminished within an elastic body due to the presence of a fracture. The fracture is modeled as having an elliptical cross section whose half-length a is much larger than its half-aperture b; the fracture cuts a plate having infinite extent in the a-b plane and unit thickness. Assuming plane strain and uniaxial extension, the strain energy released from the unit-thickness plate is (e.g., Gudmundsson, 2011) where σ is the remote applied stress from uniaxial extension, v is Poisson's ratio, and E is Young's modulus. As such, ΔU increases with a 2 . Because the fracture surface area only increases linearly with a for openingmode propagation, fracture propagation under a constant load will progressively dissipate elastic strain energy faster than the newly created surface area absorbs surface energy, and so, propagation will accelerate unstably.
Because the present fractures are actually arrays of veinlets, they generate considerably more surface area per unit increase in a. We can estimate dA/da by representing the veinlets as plane-strain fractures whose surface area, per unit thickness as in Griffith's scenario, is where L is the length of a line transecting the macrovein ellipse at a point along (−a, a) and lying at an angle θ (estimated at 70°for the calculations below) from the major axis ( Figure 8). These veinlets are spaced at a constant interval distance S along the ellipse. Equation (3) is a conservative estimation of the true surface area in that it ignores veinlet curvature in the a-b plane, but the underestimation will be small where b/S is large. As stated above, S is observed to remain roughly constant through vein growth, so the degree of underestimation should diminish with vein growth. The total A for a macroscopic calcite vein is the summation of Equation (3) across the vein: Inspection of Figure 8 reveals that the increase in the total surface area (sum of all veinlet surfaces) with a depends on the aspect ratio (b/a) taken on by the vein as it grows. We assume that this relationship is linear or sublinear, based on veins of different final sizes, as preserved in cores. For linear or sublinear scaling b ∝ a n (n ≤ 1), A increases with a 1 + n . As such, for linear b:a scaling, A increases with a 2 (Figure 9). In such a case, even if the fracture formed under a constant load, the dissipation of elastic strain energy as the fracture grew would be matched by the absorption of surface energy created by the fracture.

Journal of Geophysical Research: Solid Earth
For sublinear b:a scaling, increased surface creation during lengthening would not be sufficient to result in stable propagation under constant-load conditions but would nevertheless generate more surface area, compared to the case of an ellipsoidal crack. The best-fit power law to the observed data ( Figure 6a) indicates b scaling with a 0.75although the quality of this fit is statistically indistinguishable from that of the linearbest-fit.
We observe fractures having 2a as great as 72 mm; many veins are likely much longer but truncated by the core. At this observed length scale, and assuming linear b:a scaling, the observed veins have approximately five times as much surface area as an equivalently sized ellipsoidal fracture (Figure 9).

Implications
Natural hydraulic (fluid-pressure-driven) fractures probably do not propagate under constant loads. The source of pressurized fluid is finite, and fluid pressure probably cycled within fractures, transiently diminishing as fracture volume increased (Lacazette & Engelder, 1992). Therefore, the driving stress for natural hydraulic fractures likely diminished over time, consistent with field evidence that fractures eventually arrest. This study documents a natural example wherein propagation along the mechanical anisotropy of shales resulted in excess fracture surface area with fracture growth. This result would have produced a stabilizing effect during fracture growth, as elastic strain energy was more quickly dissipated as surface area, compared to the ideal case of two discrete fracture surfaces. We therefore conclude, from this case study of the vein distribution, geometry, and structure in the Posidonienschiefer of the LSB, that excess surface area generation can serve as a potential stabilizer of natural hydraulic fracture propagation in shales.
Artificial hydraulic fracturing must take into account fracture propagation whether as a desired result in hydrocarbon exploration or as a hazard in wastewater or CO 2 storage. Laboratory studies of fracture propagation along bedding in shales suggest that the fracture toughness, which is proportional to the square root of surface energy, may be reduced to roughly half the value across bedding (Chandler et al., 2017;Luo et al., 2018). Preexisting fractures are also planes of weakness that can deflect fracture paths (Lee et al., 2015). Sedimentary planes of weakness appear to be the reason for the stepwise propagation of the present veins. It remains to be seen under what conditions the stabilizing effect described above emerges during natural or artificial fracturing, but the necessary conditions would appear to include some mechanical anisotropy and a driving stress oriented oblique to it, promoting en échelon stepping across bedding planes.
The formation of the present fractures is attributed to fluid overpressure. Structures that are similar, in that they contain fibrous mineral fill and, in some cases, abundant host-shale inclusions (e.g., "beef" and "cone-incone"), have also been attributed to fluid overpressures (Cobbold et al., 2013;Meng et al., 2017;Zanella et al., 2015;Zhang et al., 2016). In settings not conducive to calcite precipitation, this type of fracture propagation may be common but simply not preserved. However, simultaneous calcite precipitation may have important effects on propagation. Olson (2003) suggested that hydraulic fracture propagation might be stabilized by fracture lengthening being matched with narrowing in order to preserve fracture volume, for cases in which infiltration of fluids is slow compared to fracture growth. In the case of fibrous veins, although there may be Figure 8. Representation of the surface area (per unit thickness into page) of an ideal, elliptical fracture versus that of an array of veinlets (parallel lines) that open along beddings planes and rotate the host rock to an angle θ with respect to the major axis. Variables from Equations (3) and (4) illustrated. Figure 9. Increase in surface area, per unit thickness, with length, for an ellipsoidal crack (circles) and a vein comprising a series of veinlets (triangles). Note veinlet-array fracture surface area increases with the square of length, assuming linear aperture:length scaling.

10.1029/2019JB018442
Journal of Geophysical Research: Solid Earth some small, transient fracture open space at the vein-wall interface (Hilgers et al., 2001), simultaneous cementation with fracture growth would prevent fracture closure, precluding closure as a potential stabilizing mechanism.

Conclusions
Layer-parallel calcite veins in the Posidonienschiefer of the LSB, Germany, formed as natural hydraulic fractures in response to hydrocarbon maturation-related overpressure. The veins filled with calcite cement as they lengthened and widened. Macroscopically, they lie oblique to the shale bedding fabric; however, microscopic evidence shows that they propagated by stepping across bedding planes. In fact, essentially all fracture opening occurred along bedding planes. The fractures repeatedly entrained host-rock inclusions as they grew, creating excess surface area during growth, relative to the ideal case in which fractures make only two surfaces. This excess surface area generation would have produced a stabilizing effect, by which more elastic energy is consumed per unit increase in fracture length. Thus, shale anisotropy could potentially stabilize hydraulic fracture propagation.