Mineral Fabric as a Hidden Variable in Fracture Formation in Layered Media

Two longstanding goals in subsurface science are to induce fractures with a desired geometry and to adaptively control the interstitial geometry of existing fractures in response to changing subsurface conditions. Here, we demonstrate that microscopic mineral fabric and structure interact with macroscopic strain fields to generate emergent meso-scale geometries of induced fractures. These geometries define preferential directions of flow. Using additively manufactured rock, we demonstrate that highly conductive flow paths can be formed in tensile fractures by creating corrugated surfaces. Generation, suppression and enhancement of corrugations depend on the relative orientation between mineral fabric and layering. These insights into the role of micro-scale structure on macro-scale flow provide a new method for designing subsurface strategies to maximize potential production or to inhibit flow.


.1 3D Printing
Geo-architected layered rock samples with preferred mineral fabrics were created using a ProJet CJP 360 3D printer. Layers of calcium sulfate hemi-hydrate (0.1 mm thick bassanite powders, Figure 1 left) were bonded with a proprietary water-based binder (ProJet X60 VisiJet PXL) that produced gypsum as a reaction product (Figure 1 center). The gypsum mineral growth direction is oriented by the direction of linear movement by the binder spray head that moves in a manner similar to an ink jet printer. The width of the mineral bands (see Figure 3b in the paper) is controlled by the design of the binder spray head which contains evenly spaced holes to release the binder. After one layer of bassanite is deposited, the binder is sprayed across the surface in linear rows. When a new layer of bassanite is deposited on the previous layer, the reaction of binder with bassanite powders causes gypsum crystals to form bonds between the bassanite layers. The direction of the binder spray head is an input parameter to the 3D printer. The sample geometry is designed in stereolithography CAD software (i.e. STL format) that is imported into 3D printing software. During the layout in the CAD software, the samples are oriented to produce the mineral fabric orientation shown in Figure 2 in the manuscript.

Material Properties
Powder X-ray diffraction (XRD) was performed to determine the percent bassanite (2CaSO 4 · H 2 O) and gypsum (CaSO 4 · 2H 2 O) in the 3D printed samples. The XRD system was a Panalytical Empyrean X-ray diffractometer equipped with Bragg-Brentano HD optics, a sealed tube copper X-ray source (λ = 1.54178Å), with soller slits on both the incident and receiving optics sides, and a PixCel3D Medipix detector. The anti-scatter slit (1/2 o ) and divergence slit (1/8 o ) as well as the mask (4 mm) were chosen based on sample area and starting θ angle. The XRD measurement found that the powder contained 97% bassanite while the printed samples (i.e. after application of the binder) were ≈50-50 bassanite and gypsum. The average density of the 3D printed geo-architected samples was 1190 ± 5.5kg/m 3 .  Table   1. For comparison, Table 2 provides compressional and shear wave velocities for a layered shale [1].
Unconfined compressive strength (UCS) testing was performed on 3D printed gypsum samples   Figure 3 shows the geometry of the cylinders used in the UCS testing and representative stress-strain curves for the 3D printed samples. UCS values for the 3D printed geo-architected samples are given in Table 3 and for comparison values for shale are given in Table 2. An additional comparison between 3D printed gypsum samples and natural rock can be found in Kong et al. [2]

Xray CT
A 3D X-Ray Microscope (Zeiss Xradia 510 Versa) was used to acquire 2D radiographs and to perform 3D computed tomography for the small geo-architected samples during in-situ 3PB loading test. The small samples were placed in a Deben CT5000 in-situ uniaxial loading device in the 3D X-ray microscope. The loading rate was 0.1 mm/min with load information recorded at 100 milisecond sampling rate. The settings for both the 2D and 3D scans were 80 kV and 7W Xrays,

Surface Roughness Analysis
The surface roughness maps were corrected for arbitrary rotations associated with mounting the sample in the laser profilometery system. The gradients were determined by fitting a 2D plane to the surface and then subtracting the gradients from the asperity arrays. Next, the minimum asperity height was subtracted from all points to yield asperity heights, z ( x, y) that ranged from zero to the maximum for a given surface.
The isotropy or anisotropy of a surface asperity height distribution was determined from a 2D auto-correlation analysis. In this approach, a 2D Fourier Transform, F T , of the asperity heights from a surface, Z = F T (z(x, y)), was multiplied by the complex conjugate,Z, and then an inverse Fourier transform, F T −1 , was performed on this product: and divided by the mean of the square of z(x, y) to obtain S(x, y), the 2D auto-correlation function.
The 2D asperity map was rectangular in shape which could bias or generate artifacts in S(x, y).
For each surface roughness map, the auto-correlation function, S(x, y), was calculated for 2 circular subregions (10 mm diameter). For each sample type (i.e. 3D printed samples), the presented 2D auto-correlation functions represent an average < S(x, y) > overall samples of a given type. The auto-correlation function indicates the probability that an asperity at a distance r(x, y) will have a similar height. The maximum probability is 1 when r(x, y) = 0 when a comparison is made between a point and itself.
Micro-slope angle analysis was also performed on the asperity height map, z(x, y), to serve as a measure of the relative smoothness or roughness. Park & Song [3] defined the microslope angle as the dip of the slope between neighboring asperities. A microslope analysis was performed by finding the local slope, s, where and s y = dz(x, y) dy which is the derivative of the surface roughness profile in the x-direction (horizontal and perpendicular to the direction of fracture propagation) and y-direction (vertical direction and parallel to the direction of fracture propagation). The microslope angle is taken relative to the horizontal and is found by and θ sy = arctan(s y ).
A surface was defined as relatively "smooth" if the average microslope angle distribution full-width at half the maximum,(θ save ), was θ save < 15 o . θ save ≈ 60 o has been measured on induced fractures in granite [4], while for synthetic gypsum fractures cast against sand paper, θ save ≈ 10 o , 30 o , and 40 o were observed as the grit size increased from 68 to 265 to 530 micrometers [5].
The original work by Park & Song [3] found θ save ≈ 20 o − 30 o for joints in granite and gneiss. From fracture surfaces in slate after shear box testing [6], gaussian distributions of microslope angles were observed that were centered on 5 o and −8 o for different applied normal loads with θ save ≈ 15 o , and 25 o , respectively.

Flow Simulation
The aperture distribution was obtained by placing the measured rough surface against a flat plane that creates 5% contact area to capture the effect of large scale roughness on fluid flow. This yields a N xN array of apertures that is replaced with a network or graph of elliptical pipes [7,8,9,10,11,12,13] from which a sparse system of linear equations of O(N 2 ) is generated.
For each pipe, the volumetric flow rate is given by q = ∆P/R where R is the resistance and ∆P is the head. The resistance is given by with f given by and K = a 2 h 2 . The major and minor axes of the two ellipses (which represent the ends of an elliptical pipe) that is formed between adjacent rows of apertures in the array are given by a i and b i where i = 1 and 2. Half of the maximum aperture is represented by H, while a is the average minor axis between rows. The average area of the two ellipses is given by A avg . Additional details on the approach can be found in Pyrak-Nolte & Morris [14], Cheng et al. [10], and Petrovitch [11].

Permeability Normalization
The average mean aperture and average aperture of the critical neck in the x− and y− directions are listed in Table 4 and were used to normalize the average permeability for each sample.