## ABSTRACT

Reservoir core measurements can help guide seismic monitoring of fluid-induced pressure variations in tight fractured reservoirs, including those targeted for supercritical $CO2$ injection. We have developed the first seismic-frequency “room-dry” measurements of fracture-specific shear stiffness, using artificially fractured standard granite samples with different degrees of mating, a well-mated tensile fracture from a dolomite reservoir core, as well as simple roughened polymethyl methacrylate (PMMA) surfaces. We have adapted a low-frequency (0.01–100 Hz) shear modulus and attenuation apparatus to explore the seismic signature of fractures and understand the mechanics of asperity contacts under a range of normal stress conditions. Our instrument is unique in its ability to measure at low-normal stresses (0.5–20 MPa), simulating “open” fractures in shallow or high-fluid-pressure reservoirs. The accuracy of our instrument is demonstrated by calibration and comparison with ultrasonic measurements and low-frequency direct shear measurements of intact samples from the literature. Pressure-sensitive film was used to measure real contact area of the fracture surfaces. The fractured shear modulus for most of the samples shows an exponential dependence on the real contact area. A simple numerical model, with one bonded circular asperity, predicts this behavior and matches the data for the simple PMMA surfaces. The rock surfaces reach their intact moduli at lower contact area than the model predicts, likely due to more complex geometry. Finally, we apply our results to a linear-slip interface model to estimate reflection coefficients and calculate S-wave time delays due to the lower-wave velocities through the fractured zone. We find that cross-well surveys could detect even well-mated hard-rock fractures, assuming the availability of high-repeatability acquisition systems.

## INTRODUCTION

There is increasing interest in identifying and characterizing fractures (natural and engineered) in reservoirs, particularly when fractures are the dominant form of permeability (e.g., Bear, 1972; Haimson, 1975; Zoback and Byerlee, 1975; Selvadurai and Glowacki, 2008). Most engineering efforts within fractured reservoirs (oil and natural gas extraction, enhanced geothermal systems, carbon dioxide sequestration, etc.) rely heavily on the knowledge of the reservoir permeability to conduct operations (e.g., Streltsova, 1983). Seismic imaging is the primary method used to gain detailed information about reservoir structure, which includes fracture density and orientation (e.g., Willis et al., 2006). Data obtained from field surveys must be interpreted within the framework of empirical relationships developed from laboratory tests on reservoir cores and informed by our theoretical understanding of rock elasticity and seismic wave propagation (Mavko et al., 1998). It is important to undertake a range of laboratory experiments to isolate the reservoir processes to be monitored seismically and upscale these measurements to the conditions in the field.

The bulk of laboratory seismic measurements are undertaken using ultrasonic measurement techniques developed over the past decades (e.g., Birch, 1960; Nur and Simmons, 1969; Lockner et al., 1977; Winkler and Plona, 1982; Vanorio et al., 2002; Vialle and Vanorio, 2011), by propagating elastic waves of wavelength much smaller than the sample size through the specimen to measure traveltime and thus wave speeds. This laboratory approach cannot be used at lower seismic frequencies (1–100 Hz) because the wavelengths involved are much larger than the available core specimens. Resonance techniques allow lower frequencies with long, slender samples (e.g., Birch and Bancroft, 1938), but measurements at kilohertz frequencies usually require samples with meter-scale dimensions (e.g., Winkler and Nur, 1979), unless a mass loaded resonator is used (e.g., Cooper, 1979; Nakagawa, 2011), and then they provide data only at the fundamental resonance frequency and high harmonics (e.g., Katahara et al., 1982; McCann and Sothcott, 2009). Because field measurements involve averaging over the scale of a seismic wavelength, it can be difficult to distinguish between heterogeneities at the range of scales of the field and lab and those due to intrinsic dispersion effects (e.g., Wu, 1982; Pyrake-Nolte and Nolte, 1992; Aki and Richards, 2002). A limited number of high-bandwidth, low-frequency instruments exist capable of probing rocks through the use of forced subresonance oscillation methods (e.g., Spencer, 1981; Jackson et al., 1984, 2011; Peselnick and Liu, 1987; Jackson and Paterson, 1993; Jackson, 2000; Batzle et al., 2006; Tisato and Madonna, 2012; Nakagawa, 2013; Subramaniyan et al., 2014; Tisato et al., 2014). These instruments can measure elastic moduli and attenuation in the seismic frequency range to probe processes that may cause relaxation at the lower frequencies used in the field.

Various fracture and fracture network wave-propagation phenomena, including tube waves, seismic focusing, interface waves, and guided waves, have been characterized in the laboratory and field (e.g., Hardin et al., 1987; Oliger et al., 2003; Shao and Pyrak-Nolte, 2013; Shao et al., 2015; Nakagawa et al., 2016). The size and density of fractures that can be seismically imaged in the field depend on the magnitude of the elastic contrast they induce and the wavelength of the probing seismic wave, as well as survey details such as geometry, source type, energy, and repeatability (e.g., Silver et al., 2007). The detectability of fractures can be estimated based on the reflection, conversion, and transmission of seismic waves; or from the time delay a seismic wave would experience from traveling through the slower region around a fracture. Schoenberg (1980) derives solutions for calculating the reflection coefficients from a linear-slip interface model (LSIM) dependent on fracture-specific stiffness, seismic impedance, and frequency, which have been confirmed with measurements of ultrasonic waves across fractures at varying effective stresses (Pyrake-Nolte et al., 1990). The results are frequency-dependent, but they were only initially verified in the ultrasonic range. Fracture stiffness, or compliance, has been measured at ultrasonic frequencies over a range of scales, fluid saturations, and pressures, shown to agree broadly with the LSIM theory (e.g., Trimmer et al., 1980; Worthington and Lubbe, 2007; Lubbe et al., 2008). Pyrak-Nolte and Nolte (1992) show that frequency-dependent fracture stiffness can result from the sampling of fractures with heterogeneous stiffness at a range of wavelengths, but they also point out that there could be dynamic effects, such as locking or friction. Without studying these potential dynamic effects, it may be hard to separate them from scaling effects in the field.

Nakagawa (2013) measures frequency-dependent normal and shear fracture stiffness in an unmated, fractured Berea sandstone sample (water saturated and drained) under a range of confining pressures from 0.09 to 1.7 MPa (with pore pressure under a weak vacuum) in the frequency range of 1–100 Hz. A poroelastic model was used to show that the measurements were very sensitive to partial saturation; a small amount of gas (saturation 99.5%) was needed to fit the data. To this point, there have been no fracture specific stiffness measurements of “room-dry” or mated fractures in the field seismic frequency range to determine the seismic signature of narrow single fractures in the field.

Isolating single fractures in the laboratory allows quantification of fracture surfaces and thus illuminates the importance of microscale asperity mechanics. Surface topography studies in rock have been used to predict the stress states imposed by the interaction of rough surfaces that comprise the two sides of a fracture (Brown and Scholz, 1985; Brown et al., 1986; Power and Tullis, 1991; Candela et al., 2009). Greenwood and Williamson (1966) are among the first to use statistical distributions of surface roughness to determine the heterogeneous distribution of normal stress along the fracture interface (see also Hansen et al., 2000; Persson et al., 2002; Schmittbuhl et al., 2006) using the framework of Hertzian contact mechanics (Hertz, 1882). Variations in the normal stress field are assumed to be a function of the elastic deformation of the contacting topography as they are pressed together (Johnson, 1985). Understanding the heterogeneous distribution of normal stress has proved important because it can determine local shearing strength of a fracture surface in the context of contact mechanics (Archard, 1961; Bowden and Tabor, 2001; Scholz, 2002; Campañá et al., 2011; Afferrante et al., 2012; Pohrt and Popov, 2012). A better understanding of the strength heterogeneity along fracture surfaces is crucial in predicting a fracture’s constitutive parameters (e.g., shear modulus) in natural settings.

Constitutive models determining how asperities form are complicated due to the self-affine nature of the fracture surface topography (Candela et al., 2009; Brodsky et al., 2016); individual contacts (asperities) may form due to elastic, elasto-plastic, or fully plastic interactions (Johnson, 1985) and depend on the length scales at which they are analyzed (Persson, 2006). In this study, we use pressure-sensitive film (Selvadurai and Glaser, 2015a, 2015b) to empirically measure the asperity contact formation as the fracture is subjected to increased normal stress. The pressure film allows us to measure changes in the real contact area and local variations in the normal stress field as the gross normal stress is varied.

We have adapted a low-frequency (0.01–100 Hz) shear modulus and attenuation apparatus (Bonner and Wannamaker, 1991a, 1999b) to explore the shear mechanics of fracture surfaces to help interpret field data from highly fractured reservoirs. Our instrument is unique in its ability to measure shear properties at low normal stresses, which is particularly relevant to shallow or high-fluid-pressure reservoirs. The equivalent effective stress experienced on a fracture due to fluid pressure and overburden can be calculated using a measured or inferred Biot’s (1941) coefficient. We present calibration data that show the accuracy of our instrument, as well as measurements of roughened polymethyl methacrylate (PMMA) interfaces in contact, artificially fractured granite samples (with different degrees of mating), and a fractured dolomite reservoir core (from a carbon sequestration target). The fracture surfaces are characterized with an optical profilometer and pressure-sensitive film, which measures the surface topography and normal stress distribution, respectively. The combined measurements of contact formation and their shear mechanics show an exponential dependence of the fractured shear modulus on the real contact area. We construct a simple model, which helps to explain these measurements, while neglecting the more complex effects of facture finite width, asperity geometry, and interaction (Boitnott et al., 1992; Moore and Lockner, 1995; Yoshioka, 1997; Misra and Huang, 2012; Morris, 2015).

Finally, we calculate the conditions under which these fractures should be directly detectable through seismic wave reflection and transmission. As fractures open up, decreasing modulus and $Q$ (the inverse of attenuation), with less uniaxial load, more seismic energy will be reflected and seismic velocity will decrease causing larger traveltimes. Both of these effects will increase the seismic signature in the field. In this way, we estimate the normal stress dependence of the seismic visibility of fractures (with different degrees of mating) in hard rock such as granite and dolomite. Although we are not measuring the direct effect of pore fluid, fracture measurements under room-dry conditions, in which room humidity insures some adsorbed water on grain surfaces (Clark et al., 1980), are important for understanding the stress dependence of seismic properties as well as simulating dry environments from compressed air (e.g., Majer et al., 1997) or supercritical $CO2$ injection. Dry frame properties are also important for rock-physics interpretation in general, such as for understanding the added effects of fluid substitution (Mavko et al., 1998). More measurements of natural fractures with pore fluid, incorporated into fracture network wave propagation models, are necessary to predict the seismic signature of increased pore pressure in fractures, expected with large volume injection of supercritical $CO2$ for geologic carbon storage.

## METHODS

### Shear modulus and attenuation apparatus

Our approach uses a forced torsional oscillator (e.g., Berckhemer et al., 1982; Jackson et al., 1984; Gribb and Cooper, 1998; Jackson et al., 2011), operating without wave propagation and at subresonance conditions, to directly measure the complex shear stress and strain of solid right cylindrical rock samples (9 mm diameter with various lengths; see Table ^{1}). The apparatus (Figure ^{1}) is a segmented torsional spring. One end includes a noncontact driver consisting of six electromagnets and six Nd-Fe B permanent magnets mounted on the central rotor (labeled 1 in Figure ^{1}), which twists the entire bar of the apparatus in series from a fixed point at the other end (labeled 2). Because all the components of the apparatus are assumed to be elastic on the time scales of interest, the torque (and shear stress) is transmitted evenly throughout the bar’s entire length, and by measuring the amount (approximately 10s of microradians) and timing (approximately tenths of a microsecond) of the twist at various locations along this bar, we can calculate the shear modulus and attenuation of the sample.

The magnetic driver applies a zero-frequency “pretwist” and then applies an oscillating torque around that point. This pretwist allows sinusoidal oscillating torques without twisting past the neutral point, which would add hysteresis. The oscillator operates at frequencies between 0.01 and approximately 100 Hz, with the upper limit dictated by system resonances that depend on the sample length and stiffness. The driver is capable of achieving strains ranging from $10\u22127$ to $10\u22124$, spanning the linear and nonlinear domains of the specimen.

The fixed end (labeled 2) of the apparatus is attached to a hollow section of aluminum (labeled 3), with arms and targets at each end (labeled 4) for eddy current proximity detectors (labeled 6; Kaman Sciences KD 2300-0.5SU and KD2300-2S). The arms on both sides of the specimen allow for measurement of specimen twist. The mounted arms (labeled 4 and 5) are made of a carbon fiber honeycomb material, to minimize weight and maximize stiffness. The detectors (labeled 6) are configured to measure differential displacements, which depend only on the elastic deformation of the known aluminum segment (labeled 3). These sensors (labeled 6) are used to determine the shear stress being applied throughout the apparatus by the magnetic driver (labeled 1). A set of aluminum collets (labeled 7; Acura-Flex collet chucks with straight shanks) holds the rock sample (labeled 8), attaching it in series with the aluminum bars (labeled 3) of the apparatus on each side. The twist is again measured on the driver end (labeled 5) of the sample, which captures the shear strain in the rock.

By taking a Fourier transform of the measured signals using a dynamic signal analyzer (Stanford Research Systems SR785), we can isolate the response at the driving frequency and obtain the complex shear stress and strain. Following classical anelasticity theory (Zener, 1948; Norwick and Berry, 1961; Gordon and Davis, 1968; Lakes, 2004), the ratio between real parts of the stress and strain gives the shear modulus and the tangent of the difference in phase angles (imaginary parts) provides the shear attenuation. We process the stress-strain curves for a range of amplitudes and frequencies to calculate the shear modulus and analyze the linearity and dispersion properties of the rock samples. Signal averaging is used to improve the signal-to-noise ratio; this is especially important for measuring phase angles accurately. All data acquisition and control functions are automated using a National Instruments LabView code, written for this instrument.

In torsion, shear strains are not uniform throughout the cylinder of rock, and increase linearly with radius. The outside of the sample experiences the largest strains, thus it is overrepresented in the measurements of macroscopic modulus and attenuation. For this reason, we calculate and present the maximum strain values experienced at the outside surface of the sample. As long as these maximum strains are lower than approximately $10\u22126$, the cross-sectional variation in strain in the sample remained in the linear regime (verified by the linear stress-strain relationship) and should not bias our results. The cylindrical geometry is also taken into account when calculating the shear stress by using the moment of inertia for a solid cylinder (see equation ^{1}).

Unlubricated hard bearings inside the driver (labeled 10) suppress flexural modes, ensuring we deform the samples purely in torsion. There is also a uniaxial loader (labeled 11), a spring-loaded fixture with another hardened ball bearing (labeled 10) at the far end of the driver that allows the application of up to 20 MPa of uniaxial stress along the centerline of the apparatus. Bellville washers and a compression nut (labeled 12) generate the restoring force, ensuring that it is linear with normal displacement. A load cell (labeled 13), including a calibrated strain gauge, measures the axial strain of a known aluminum piece, and it is calibrated to provide the applied uniaxial stress. This loader provides a normal stress to close cracks oriented in the plane of the cylinder’s diameter (those that most affect the torsional response). The hardened bearing (labeled 10) provides a point contact to minimize the torsional damping due to the uniaxial loading assembly, making it small compared with that of a fractured sample; use of this feature of the apparatus still raises the lower bound of attenuation we can reliably measure, which will be addressed later in the “Error analysis” section (Appendix ^{2}).

The instrument is well-suited for the study of fractured rocks for several reasons. The uniaxial stress assembly provides normal stress to close up fractures and slowly open them with decreasing normal load. This can simulate, using Biot’s (1941) coefficient, important processes in fractured reservoirs, such as effective stress changes from injection or production. However, our study does not directly probe matrix/fluid interactions, such as Biot’s (1956) or squirt-type (Mavko and Jizba, 1991) effects. Our apparatus is particularly well-suited for very low normal stress measurements, which is inherently difficult for other instruments that rely on confining pressure of at least 10 MPa, and jacketing, to make frictional contact with the sample (e.g., Jackson et al., 2011; Li et al., 2014).

### Calibration and error estimation

The stress and strain measurements are made by eddy current proximity sensors, calibrated using known values of torque and displacement. The calibration of attenuation is more straightforward because it is independent of the sample dimensions. Details of the calibration techniques and equations for calculating shear stress and strain from our raw measurements are given in Appendix ^{1}.

There are two different sources of error in our measurements: random error in the measurements themselves due to electrical and numerical noise and systematic errors from imperfect calibration. We minimize the electrical and numerical noise by taking 1000 repeated measurements at each stress condition with all of the samples. Each source of systematic error has also been quantified and compared with the variance in the measurements; both types of error are described in Appendix ^{2}. We choose to give error bars based on the standard deviation because this gives a better comparison between our measurements of the same rock sample under different stress conditions.

### Surface characterization

To better understand the micromechanics of our fracture surfaces and how they change with normal stress, we used a novel stress distribution imaging technique (Selvadurai and Glaser, 2015b). A Fujifilm Prescale medium range (12–50 MPa) pressure-sensitive film measured the contacting asperities formed between the interacting surfaces under a range of applied normal stresses. Details of the film’s capabilities and calibration are given in Appendix ^{3}.

For geometric comparison, surface topography was also directly measured using an optical scanning profilometer (Nanovea $PS50/3000\u2009\u2009\mu m$ optical pen). The spatial resolution in the $r\u2010\theta $ plane was $50\u2009\u2009\mu m$ along the fracture surface and $0.5\u2009\u2009\mu m$ in height. The surface topography of the fractures varies significantly between our samples, depending on grain size and fracturing method, and the scans are useful for characterizing the tilt of the fracture and surface roughness. A profile across the center of the sample diameter is used to describe the topographic trend of the surface (either tilted, bowed, or with multiple ridges/valleys), as well as a sense of the surface roughness (Figure ^{2}). These profilometer images are not in the same orientation as the film images shown below. Although the initial contact conditions measured from the pressure-sensitive film should ultimately determine the physical properties of the fracture (e.g., shear stiffness; Berthoude and Baumberger, 1999), the surface topographies were measured to provide a better understanding of the fracture and help to validate some of the assumptions made when using the pressure-sensitive film.

## RESULTS

### Comparison of standards

We compare our measured modulus values with those obtained by the ultrasonic pulse transmission method for intact samples of aluminum, PMMA (or acrylic), ABS plastic, Duperow dolomite, and Westerly, Sierra White (SW), and Montello granites (Table ^{1}). Ultrasonic tests were conducted using 5 MHz S-wave transducers (Panametrics V155), minimal manual contact pressure and ultrasonic couplant, and subtracting the “face-to-face” piezoelectric rise time. The ultrasonic waves were shot across the diameter of the cylindrical samples, the transducers inline contact with the sample, to avoid TI anisotropy in aluminum from the extrusion process, which the torsional measurements average. The moduli measured at low frequencies are systematically lower than at ultrasonic frequency for the more attenuating samples, which also show some dispersion within the frequency range of our apparatus. PMMA has well-known dispersion due to unwinding of complex polymer chains (Ferry, 1980), whereas room-dry rocks and metals can show dispersion from adsorbed water or other viscous/frictional regions on microcracks or grain surfaces (Zener, 1948; Cooper, 1979). In contrast, aluminum 6061-T6 shows no difference within the error bars between measurements at different frequencies because it has negligible attenuation as discussed in Appendix ^{1}. The difference in modulus between measurement techniques is likely due to a combination of these dispersion characteristics of each material and the added importance of the outside of the sample in our measurements, due to the torsional configuration’s increased peripheral strains determined with equation ^{2} (in Appendix ^{1}). The outside rind of the sample is likely to have a lower modulus than the sample in bulk, because of open microcracks generated during the coring process; these cracks are likely to be open because all of these comparison measurements were made without uniaxial load. Our room-dry, atmospheric pressure, low-frequency shear-modulus measurements of Westerly granite (WG), 17–19 GPa, are consistent with the scarce available literature data of comparable conditions. Cooper (1979) measures the torsional resonance at frequencies between 0.649 and 0.685 Hz, giving shear moduli of 16–18 GPa, in different directions to show its slight anisotropy. Simmons and Brace (1965) find even larger (approximately 50%) discrepancies beyond the expected error between dynamic and static measurements at atmospheric pressure, attributing them to open cracks that effect strain-gauge measurements across the entire sample more than ultrasonic wave propagation, which may bypass cracks.

We found that PMMA provides the best low-frequency standard when compared with published modulus and attenuation results from a range of historic and modern systems (Koppelmann, 1958; Yee and Takemori, 1982; Capodagli and Lakes, 2008; Nakagawa, 2011; Tisato and Madonna, 2012). Although this material often varies in absolute mechanical properties from batch to batch due to variable fabrication procedures, the overall frequency dispersion characteristics of the material are fairly reproducible between batches. In Figure ^{3}, we compare our measurements of PMMA with those published by other low-frequency or wide-bandwidth instruments. Noting that absolute shear moduli of PMMA can vary significantly between batches, the similarities in the dispersion characteristics we measure with those measured in the literature provide confirmation of our system and calibration protocol.

### Fracture results

In this section, we document some of the first results showing the stress-dependent moduli of fractured materials at seismic frequencies and normal stress states between 0 and 16 MPa. The measurements were made at the highest normal stresses first (often approaching intact behavior), with subsequent lower stresses representing the fracture opening.

Besides serving as a standard for comparing our measurements with literature values, we also used PMMA to test our ability to probe the elastic properties of fracture surfaces. To simulate a roughened fracture, the ends of our PMMA samples were sanded with P80 grit emery cloth (approximately $201\u2009\u2009\mu m$ particle size). A single piece of PMMA was used to measure the intact modulus and attenuation, shown by blue circles in Figure ^{4} and ^{4}, respectively. Next, two pieces of PMMA were combined to study the roughened “fracture” surface between them. This measurement gives the fractured shear modulus and attenuation, shown by red circles also in Figure ^{4} and ^{4}, respectively. By subtracting the shear compliance, $\eta =1/G$, of the intact sample from the shear compliance of the fractured sample at every level of uniaxial stress, using the principle of superposition for linear elasticity, we can estimate the shear compliance of the fracture itself (Jaeger and Cook, 1969). In Figure ^{4}, we show measures of PMMA shear compliance for the intact sample (blue circles), fractured sample (red circles), and of the fracture only (green circles). Because PMMA is relatively compliant and attenuating compared with rock, the effect of the fracture is somewhat masked. The fractured sample reaches the modulus of an intact sample at high normal loads, implying a negligible compliance for the fracture itself (Figure ^{4}). As the fracture opens due to decreasing normal stress, there is a clear effect on shear modulus (Figure ^{4}) and shear attenuation (Figure ^{4}). The frequency dispersion does not change significantly with normal stress; therefore, we only graph the value at 8 Hz for each stress condition.

The closing and stiffening of the fracture with normal load correlate with the increase in real contact area ($Ar$), as calculated from the pressure film (Figure ^{5}). Local normal stress distribution was also measured with the film; the images for uniaxial loads from 1 to 9.5 MPa are shown in Figure ^{5}–^{5}. The contact surface was dominated by the slightly beveled topography of one of the surfaces, giving a growing circular contact at the center that grows linearly with normal stress.

To investigate the behavior of a mated rock surface, we measured the properties of a WG sample before and after inducing a fracture. We first measured the intact shear modulus and attenuation of the sample, while compliance was calculated, each plotted in blue in Figure ^{6}–^{6}. No measurements were taken for stresses between 0 and 4 MPa while the specimen was intact because it was uncertain whether the fractured sample could be tested at such low normal stress without the fracture slipping. Low-stress values for the intact sample were estimated with linear extrapolation from the higher stress measurements trend. We introduced a fracture perpendicular to the centerline of the cylindrical sample at a central location of the active length (i.e., $L/2$). The periphery of the sample was scored at a central location, and then a fracture was propagated using three-point bending with a hydraulic jack assembly. Although it is commonly assumed that a leading zone of microcracks initiates fracture nucleation in brittle rocks (Reches and Lockner, 1994), almost no material was lost from our fracture faces during the fracturing process, suggesting that the damage zone of precursor cracks is on the order of a grain diameter in thickness. By using the same sample for intact and fractured measurements, the calculation to remove the compliance of the intact segments of the sample becomes more accurate. At lower normal stresses, the intact granite is considerably stiffer than the fracture, with a higher $Q$ (Figure ^{6} and ^{6}). The granite sample also shows little change in dispersion over the range of frequencies; the reported value was measured at 8 Hz.

The topography of the WG surface was much rougher than the PMMA sample, due to the grain-scale heterogeneity of the rock and the fracturing approach, as can be seen in the optical topography scan (Figure ^{2}). Figure ^{7} shows the results from the pressure-sensitive film used to characterize the WG fracture. The results are presented in a manner similar to Figure ^{5} that described the PMMA-PMMA interface. Figure ^{7} displays the relationship between the normal stress and the measured ratio of real contact area to nominal contact area ($Ar/A0$). The inset images (Figure ^{7}–^{7}) show processed pressure film scans that display the normal stress distributions (assuming a planar fracture) for various uniaxial stresses mentioned in the caption. The increased topography is mostly visible as a ridge along the bottom of the sample, which dominates the contact area throughout the range of normal loads measured (Figure ^{7}–^{7}). This ridge provided interlocking surfaces that are not perfectly perpendicular to the cylinder, possibly translating some of the shear stress into normal stress and helping to hold the fracture together. Our induced fracture was very well-mated and reached completely clamped conditions, same shear modulus as for the intact sample, at normal stresses of 7.5 MPa (Figure ^{6}).

In contrast, we measured an existing fracture in SW granite. This granite is coarser grained (approximately 1 versus 0.2 mm), and the fracture was not as well-mated or interlocking as the induced fracture in WG. Figure ^{8} shows the intact values of the shear modulus (blue circles) and the fractured modulus (red circles). We see that the fracture shear modulus never recovers to the intact modulus in the experimental range of uniaxial normal stresses. Because the fracture now contributes a measurable compliance even at the highest normal loads, we now show the calculated fracture-only shear modulus in green in Figure ^{8}. In contrast, the shear attenuation for the fractured SW granite does reach the intact value at high normal stress, and it increases gradually with lower normal stress until it increases dramatically at approximately 4 MPa (Figure ^{8}).

The relationship between real contact area and normal stress measured from the pressure film for SW granite is shown in Figure ^{9}, and visual distributions of local normal stresses are shown in Figure ^{9}–^{9}. Unlike the other samples, the SW sample preferentially formed new real contact area at the center of the interface (Figure ^{9}–^{9}), where the shear strains are lower in our testing configuration as described by equation ^{2}, and they thus have less impact on the fracture’s shear modulus.

As an example of a carbonate reservoir core, we measured an artificially fractured dolomite core from the Duperow Formation, a carbon sequestration target in north-central Montana. The sample was cored from the Danielson well (API 811151), at a depth of 1017.3 m (3337.7 ft). An adjacent plug (34B) was measured to have a density of $2.716\u2009\u2009g/cm3$, permeability of 0.01 mD, and porosity of $2.51\xb10.07%$ (Kazimierz et al., 2004; TerraTek, 2014). We also measured the ultrasonic P- and S-wave velocities for the plug at $5340\xb1160$ and $3300\xb1100\u2009\u2009m/s$, respectively. The cores, as well as preliminary surface seismic and well-log data, suggest that this reservoir has relatively low matrix permeability but many healed and partially healed natural fractures. Preliminary modeling (Zhou et al., 2013) suggests that the pressure change will be significant and potentially the largest seismic effect from the proposed $CO2$ injection; existing open fractures should be sensitive to such increases in pore pressure. By studying a tensile fracture in the reservoir core, we estimate the effect changes in stress will have on seismic reflection and velocity from similar fractures in the field.

The fracturing technique for the Duperow dolomite differed from the WG. A custom-machined holder created the tensile fracture, in which 12 sharpened screws are positioned into a groove carved around the diameter of the sample. As the screws were slowly tightened evenly around the diameter, they forced open a tensile fracture that remained perpendicular to the diameter of the cylindrical sample. The result was a well-mated, fine-grained tensile fracture that had a smoother topography than the granite fractures created by bending (Figure ^{2}). Studying the shear properties of the dolomite was done using the same approach as the previous samples; the intact modulus and attenuation were measured at various uniaxial loads (Figure ^{10} and ^{10}). The dolomite’s intact, fractured, and fracture-only shear compliances are presented in Figure ^{10}.

Figure ^{11} shows the pressure film measurements for the Duperow dolomite fracture. The relationship between the real contact area and uniaxial stresses is shown in Figure ^{11}, whereas the insets (Figure ^{11}–^{11}) show snapshots of the local normal stress distribution at the indicated normal stresses. The relatively flat surface contributes many evenly spaced asperities across the sample cross section, which converge into larger asperities at higher normal stresses. Based on our shear measurements, we find that the fracture was almost entirely closed until an approximately 6 MPa normal load was applied, but it opened gradually, lowering the shear modulus and $Q$, with decreasing normal stress (Figure ^{10} and ^{10}). This behavior suggests that the fracture was well-mated, similar to the WG fractured sample discussed previously, even though it had less interlocking topography.

## DISCUSSION

### Numerical model

Most of the fracture results (except from the SW granite) follow a consistent pattern, suggesting that the relationship between real contact area and shear modulus could be described with a simple model. We attempted to model this relationship with the simplest configuration possible: a sample with a real contact area on the fracture represented by an infinitely thin circular bonded area in the center of the cross-sectional area.

Figure ^{12} shows the general configuration of the model, and the numerical mesh is shown in Figure ^{12}. The boundary conditions and numerical schemes used to simulate the application of twist to the ends of the cylindrical sample are described in Appendix ^{4}. Using the numerical model, we calculated the torque $T$ at the ends ($z=0$ and $L$) generated from the application of the boundary conditions described by equations ^{1}–^{7}. The error in the numerically calculated torque $T$, which can be used to estimate the maximum theoretically shear stress $\tau max$ (at $r=r0$) using equation ^{1}, was only $\u22120.05$ to 0.2% of the maximum shear stress calculated numerically.

^{12}, the red stars represent the numerical results using material properties measured for the PMMA (Table

^{1}): $L=58\u2009\u2009mm$, Poisson’s ratio $\nu =0.32$, and intact shear modulus of $G=1.64\u2009\u2009GPa$, each of which were measured independently. The numerically calculated torque $T$ was determined for various ratios of $Ar/A0$ and converted to a fractured shear modulus using equation

^{3}. We see that as the ratio of real to nominal contact approaches one (i.e., $r*\u2192r0$), the model converges to the intact shear modulus $G$ and, conversely, as $r*\u21920$, the shear modulus tends toward zero. For this reason, we attempt to characterize the numerical model using an exponential relation of the form

^{4},

^{6},

^{8}, and

^{10}) and $\lambda $ is a variable that describes the manner in which stresses are transferred across the discontinuity (interface) and is likely determined by the sample geometry and the Poisson’s ratio in the model. In the experiments, fractures with more heterogeneous, out-of-plane, or interacting asperities reach higher moduli at lower real contact areas, which is captured by a higher $\lambda $ fitting parameter in this formulation.

In Figure ^{12}, we show that the numerical model (red stars) for the PMMA parameters could be fitted using equation ^{2} (black line) using a value of $\lambda =5.596$ ($R2=0.99$) and assuming $G=1.64\u2009\u2009GPa$. Comparing these results with the experimental data from the PMMA tests (blue circles), we see that the model provides a reasonable understanding of the evolution of the fractured PMMA shear modulus for a range of real contact area. The effectiveness of the model may be due to the rounded high in the center of the sample’s surface topography (Figure ^{2}) causing the circular real contact area to increase radially outward as the uniaxial normal stress is increased (Figure ^{5}) — geometrically similar to the numerical model described in Figure ^{12}.

To determine if the model could successfully replicate rock-rock interfaces, we attempted to model the dolomite experiment in a similar manner. A new numerical model was formulated to match the active length of the dolomite samples $L=25.8\u2009\u2009mm$. A Poisson’s ratio of $\nu =0.1$ (measured at ultrasonic frequencies) and the intact shear modulus of $G=16.13\u2009\u2009GPa$ (measured under the highest normal stresses) were prescribed. The relationship calculated for the dolomite parameters is different from that for PMMA (red stars in Figure ^{12} versus 12c): The dolomite model requires higher levels of real contact to recover the intact shear modulus (at $r*=r0$), characterized by a lower value of $\lambda =3.651$ ($R2=0.9804$). The calculated relationship also differs from the experimental results from the Duperow dolomite (magenta circles) characterized by $G=16.13\u2009\u2009GPa$ and $\lambda =14.13$ ($R2=0.9688$). The comparison of numerical with experimental results shown in Figure ^{12} indicates that the dolomite shear modulus approached its intact modulus at lower real contact areas than the numerical model predicted.

In Figure ^{13}, we plot the relationship between fractured shear modulus and fractional real contact area for PMMA, as well as the well-mated WG and Duperow dolomite. The exponential behaviors of the WG and Duperow dolomite were quite similar, and neither would be well-characterized by the model that was able to predict the behavior of the PMMA surface. The WG and Duperow dolomite increase modulus more with increasing contact area than the model predicts; this is likely due to the effects of the heterogeneous rock surface, finite-width damage zone, asperity interaction, and nonnormal, “shear interlocking,” fracture surfaces (Selvadurai and Yu, 2005). Surface topography has been shown to affect the constitutive behavior of fractures, including a variety of behavior due to the ratio of normal to shear fracture stiffness and local asperity interactions (Boitnott et al., 1992; Selvadurai and Yu, 2005; Campañá et al., 2011; Mirsa and Huang, 2012). Local stress transferred through these asperity interactions could cause the fracture surface to behave as if it has a higher contact area. These effects could be accounted for with a model that includes the complex fracture surface geometry (e.g., Morris, 2015), which could be input from the measured topography (Figure ^{2}) or images of stress distribution at increasing normal load (Figures ^{7}, ^{9}, and ^{11}).

The model assumes that the fracturing process represents the creation of a single, infinitely thin, through-going fracture between the two sides of intact rock. This ignores the effect of the process zone of precursor microcracks, some of which coalesced into the fracture, whereas others were left isolated on either side of the fracture surface (Reches and Lockner, 1994). As stated in the fractured results section above, the fact that few grains were broken off from fracturing suggests that these microcracks extended on the order of a full grain diameter (approximately 0.2 mm for WG) into the rock surfaces. Moore and Lockner (1995) observe extensive microcracking approximately 10 mm away from a controlled, not through-going, shear fracture in WG, under much higher stress conditions: 50 MPa of confining pressure. Direct observation or modeling of this damage zone is beyond the scope of this work, but we believe that the microcracking was much less extensive in our case because the tensile fracture was under much less stress, limited only by the tensile strength of the rock, and because it propagated much faster through our entire sample. Still, the existence of these microcracks in the rock fracture zones (as opposed to PMMA) could also help to explain why the rock fractures cannot fit the modeled contact area versus modulus relationships (Figure ^{12}). The fractured rock modulus will increase more with normal stress as these other cracks close up, adding to the effect from increasing the contact area that we model. For reference, the modulus change due to stress-induced microcrack closure is shown in the stress dependence of the intact rock sample (blue lines in Figures ^{6}, ^{8}, and ^{10}). This represents a lower bound on the magnitude of the effect, as the microcrack density has been found to increase by a factor of three or more in the damage zone (Moore and Lockner, 1995).

The topography of the WG surfaces and Duperow dolomite to a lesser degree were significantly rougher than the PMMA surfaces, as seen quantitatively by the height variations shown in Figure ^{2}. We believe that was due to the grain-scale heterogeneity of the rock and the different fracturing approaches. The induced fractures appeared very well-mated and reached completely clamped conditions (i.e., same shear modulus as the intact sample), at normal stresses greater than or equal to approximately 7.5 MPa for the WG and 6 MPa for Duperow dolomite. In contrast, an existing fracture in SW granite was not as well-mated or interlocking as the induced fractures. The SW was the only fractured sample that was limited to a shear modulus of approximately 50% of the intact shear modulus at the highest normal stresses (Figure ^{8}). It was also the only sample that seemed to preferentially add contact area close to the center of the sample cross section (Figure ^{9}), in which the shear strains are lowest and thus affects the fracture shear modulus the least. Our model could be adapted to match the radius at which contact area is added, as measured with the pressure-sensitive film, which should affect the rate at which the fractured shear modulus starts to asymptote toward the intact value. The SW behavior also does not fit an exponential relationship because it appears to show different stiffening regimes at different normal loads. Less well-mated rock fractures likely need a more advanced model to account for frictional partial slip (Cattaneo, 1938; Ciaveralla, 1998; Selvadurai and Glaser, 2015a); they do not satisfy the assumption of completely bonded contact in the low normal stress regime.

### Detectability calculation

Following Pyrake-Nolte et al. (1990), we calculated the reflection and transmission coefficients for an S-wave normal incident on our dry fractures over a range of normal stresses at seismic frequencies. They derived expressions for these coefficients using an LSIM (Schoenberg, 1980) with a single, infinitely thin fracture represented as a displacement discontinuity between two elastic half-spaces (of seismic impedance $Z$). The model does not include the effect of other parallel microcracks, which would add to the reflectivity of the damage zone. It also assumes that stress is continuous across the fracture but that displacement is discontinuous. This is a reasonable assumption for our experiment for all but the most open fractures (lowest normal stress conditions on less-mated SW), where we observe partial slip and the stress transmission is incomplete and unpredictable across the fracture. At this point, the fractures are very compliant and should reflect enough energy to be detectable. Analysis of the highest normal stress measurements was also problematic; the calculation of fracture specific shear stiffness is inaccurate when the fractured modulus is close to the same as the intact modulus because the fracture approaches infinite stiffness. These fractures have very small apertures, and reflect little seismic energy, thus they are not the focus of this study.

We calculated fracture specific stiffness values $\kappa \theta $ from our measurements by subtracting out the compliance of the intact sample, the same calculation described above. This allowed us to effectively estimate the displacement on the fracture itself for a given applied shear stress. We believe this is equivalent to a measurement of stress-displacement curves commonly used to calculate fracture specific stiffness (Pyrake-Nolte et al., 1990). To convert our specific stiffness values to commonly used units ($Pa/m$), we divided by the circumference of the sample, which converts our shear strain back to displacement on the outside of the sample radius; strain is based on an angular displacement (see equation ^{2}) only across the fracture. We used shear modulus measurements of the intact rock to calculate the impedance of the elastic half-space in the seismic frequency range.

Pyrake-Nolte et al. (1990) also derive an extension of the model to capture fluid-filled fractures with viscous coupling and dissipation, in which there is a discontinuity in velocity as well as displacement across the fracture. Although we do not directly address this scenario here, these expressions can be used to see the effect of different fluids using our fracture measurements, by assuming that the fluid itself has no effect on the fracture stiffness in shear.

Although we are interested in the seismic signature of many fractures in a reservoir as effective stress changes with injection, we analyze the effect of a single fracture from our measurements. Further analysis is necessary to understand the effect of a natural or engineered fracture network and other details of wave propagation in this setting. Figure ^{14} shows the calculation of reflectivity as a function of normal stress for each of the fractured rock samples, with the SW fracture shown at each frequency measured. The SW fracture becomes more reflective at higher normal stresses and reaches a larger value of reflectivity than the other materials, likely because it is less well-mated and thus the joint is more compliant.

Another way to determine the effect of fractures on seismic measurements is the direct time delay of S-wave arrivals due to the lower velocity of the fracture and damage zone. The time delay was calculated using measured modulus and active length to give a traveltime across the entire fractured sample, including any damaged zone, minus that in the intact sample. Figure ^{14} presents the time-delay calculation for each rock sample (there is very little frequency dependence) as a function of normal stress. Detectability of time delays as low as 6 ns has been achieved with well-calibrated semipermanent deployments of repeatable sources at higher frequencies in a crosswell configuration (Silver et al., 2007) or at low frequencies from surface sources (Dou et al., 2016). Even the well-mated single fractures at our highest normal stresses have calculated time delays more than two orders of magnitude above this instrumental limit. In practice, most deployments are not capable of such sensitivity but microsecond repeatability is straightforward to measure with piezoelectric or magnetostrictive sources. These calculations are only meant as an illustration for how these measurements could be used to guide field predictions and interpretation; more detailed wave propagation models through fracture networks and further measurements of natural fractures are necessary to more accurately predict the seismic signature in a fractured reservoir.

## CONCLUSION

We present a technical description of a low-frequency (0.01–100 Hz) shear modulus and attenuation instrument adapted to measure fractured rock samples in a torsional configuration. We provide details of the calibration and error estimation methods using a sample of PMMA that conforms to shear measurements made in the literature. Our apparatus is well-suited for the measurement of fractured samples under a range of normal stresses, from fully clamped (intact) to nearly open fractures. We observe very different behaviors for two granite samples that have different surface topography and degrees of mating as measured with an optical profilometer. We highlight the importance of the out-of-plane topographic geometry in the shear behavior of a fracture by comparing the experimental results with a simple numerical model. We combine our measurement of elastic changes due to normal load with a novel surface characterization technique that uses a pressure-sensitive film to image the normal stress distribution on the fracture interface. This provides a better understanding of the effect of microstructure and heterogeneity on seismic properties.

A simple numerical model is able to explain the exponential relationship between the real contact area and the fractured shear modulus. We were able to match the data from PMMA interfaces with the model because of their relatively simple surface topography, whereas the rock specimens exhibited more-complex behaviors. The well-mated fractures in WG and Duperow dolomite specimens followed an exponential relationship but stiffen more at a lower contact area than the model suggested. We propose that the out-of-plane asperity geometry and interaction for the more-heterogeneous surfaces become important when understanding the shear and seismic behavior of fractures at small length scales.

We then apply our fracture measurements to theory developed to calculated reflection and transmission coefficients for normally incident seismic waves in a field setting. These calculations, as well as estimations of S-wave time delays due to the fractures, are useful for determining under what stress and fracture mating conditions a fracture is seismically visible. These results are limited, representing only the specific, single fractures we measured, which are likely better mated than corresponding fracture networks in the field. Further measurements of natural fractures from the reservoir and more detailed wave propagation models through realistic fracture networks can help guide seismic monitoring, such as needed during injection of supercritical $CO2$ for geologic carbon storage.

The images of stress distribution on fracture faces in combination with sensitive measurements of large modulus changes at low normal loads provide a detailed look into underlying mechanics of seismic waves incident on fractures. Further experiments are needed to probe the range of phenomena present in fracture systems, especially at low normal loads. These laboratory measurements could be used to extract more information from open fractures in the field, but understanding the results will require further model development, both capturing the geometry of asperity contacts and frictional partial slip on contacts at lower normal stresses.

## ACKNOWLEDGMENTS

This research was funded as part of the Big Sky Carbon Sequestration Partnership (BSCSP) by the U.S. Department of Energy and the National Energy Technology Laboratory through award number DE-FC26-05NT42587. We would like to thank L. Spangler and the BSCSP leadership for access to the Duperow sample discussed in this paper. We would also like to thank S. Nakagawa, who provided useful discussion and essential characterization help in the laboratory, R. Lakes, who provided access to data, B. Buffet, who provided useful suggestions and discussion, D. Swantek, who contributed to the instrument schematic, C. O. Boro, who provided design and fabrication expertise during the early stages of the effort, and I. Jackson and an anonymous reviewer, who provided valuable reviews. P. A. Selvadurai would like to acknowledge the support of the National Science Foundation (grant no. CMMI-1131582) and the National Science and Engineering Research Council of Canada (grant no. PGSD3-391943-2010).

### CALIBRATION

The same method was used to calibrate the applied uniaxial load measurement. A purpose-built aluminum load cell was loaded in series with the apparatus, and strain was measured with four electrical resistance strain gauges in a bridge configuration (Micro-Measurements EA-06-125TQ-350). The strain measurement was calibrated to determine a linear relationship between the uniaxial load and strain gauge output. The uniaxial load was divided by the sample cross section that it was applied over to give a stress value (MPa). The linear relationship ($R2=0.99$) over a range of 30 MPa was obtained by measuring the strain with a series of known weights as well as a hydraulic press with a range of 227 kg as determined with a calibrated Bourdon tube pressure gage.

Attenuation was determined by the phase lag between torque and twist (Gordon and Davis, 1968), which is independent of sample size. Using the small-angle approximation, we directly took the phase delay (in radians), instead of the tangent of the phase delay, between stress and strain signals at the oscillating frequency as the attenuation in the sample. There is attenuation in the instrument (system $Q$), but this represents the minimum attenuation measurable, whereas most attenuating rock samples are well above these values, particularly at low normal stress states. We obtained an estimate of this minimum attenuation by measuring a low-attenuation sample, aluminum alloy 6061-T6 (for reference, $Q$ of another aluminum alloy 24ST was measured as $2.5\xd7105$ at 840 Hz by Zemanek and Rudnick, 1961). We found that the instrument attenuation was higher at low frequencies ($<2\u2009\u2009Hz$) because of low-frequency electronic drift and $1/f$ noise, as well as at low strains ($<10\u22126$), where the signal-to-noise ratio was poorer. Even for those measurements, the system phase delay was less than 0.001 radians ($Q$ equals 1000). This is at least an order of magnitude less than rock attenuation values for samples of interest, particularly for fractured rocks.

### ERROR ANALYSIS

Repeated measurements of modulus and attenuation follow a normal distribution, as evaluated by the Anderson Darling test (Anderson and Darling, 1954; see Figure ^{1} for histograms of repeated modulus and attenuation measurements). Modulus has a significantly smaller standard deviation than attenuation, as a percentage of the value. The standard error on the median of these distributions decreases with the square root of the number of samples taken, improving the reproducibility of our measurements by taking 1000 samples at each frequency; this is especially important for the attenuation measurement. This number was chosen because standard deviation appears to asymptote at approximately 1000 measurements. These measurements cover a range of amplitudes and aid in assessing the linearity of the stress-strain curves.

The linear least-squares method gives an error on the slope for each line in our calibration; by propagating this error, we can see the relative effect of each source on the measurement of shear modulus made using equation ^{3} (from Appendix ^{1}).

The hanging weight calibration of shear stress was performed nine times to improve the statistics and contributes a $\xb10.104%$ error to the shear modulus. The strain sensors were calibrated with a purpose-made micrometer, which moves the eddy current sensor in known increments relative to a metal target, and contributed $\xb10.062%$ to the modulus error. Both of these errors (together $\xb10.166%$) can be improved by making repeated measurements, but the error due to measurement of the active length $L$ of the calibration sample is more difficult to minimize. The sample was held by tightened aluminum collets, and the active length $L$ was the part of the sample that is actively deforming, which can differ slightly each time the sample is loaded due differences in how far it sat in the collets. To minimize this source of error, the sample was only unloaded and reloaded when necessary, such as for fracturing. We measured the active length by marking where the sample fits in the collets while it was loaded, and then measuring the deforming length of the sample between the collets after the sample was removed. There is some added uncertainty in the active length due to edge effects; it has been suggested that the effective length is less than the distance between clamping points (Cooper, 1979). Given these difficulties, we estimate 93.4% accuracy in measuring the “active length” of the smallest calibration samples (error of $\xb10.127\u2009\u2009cm$ from a length of 1.93 cm), which propagates to 6.6% error in the calculation of modulus for these smallest samples (the error is much smaller as a percentage for the longer samples). The measurement of the active length is also a part of the zero-length offset calculation, but it only contributes an approximately 3% error to the modulus through this calibration.

Another source of systematic error is the collets. We are careful to tighten the collets to the same torque with every loading. If they are not sufficiently tightened, the sample will slip within the collet, which is immediately apparent in the amount of strain recorded. If the collets are tightened too much, they will potentially open microcracks, or create a through-going fracture in the rock sample, also apparent in the measurements. For all of these reasons, it is crucial that the loading is done very carefully and as infrequently as possible. There are further environmental sources of variability, such as humidity and temperature changes, which are only partially controlled in our laboratory environment, and they are included in the variability of the 1000 measurements greater than 1.5 h that we average.

### PRESSURE-SENSITIVE FILM

The pressure-sensitive film is polyethylene-based and has a thickness of approximately $90\u2009\u2009\mu m$. The film has embedded microcapsules, which have a spatial resolution and were digitized to $20\xd720\u2009\u2009\mu m$. When compressed, the capsules crush and release ink with colors proportional to the applied pressure (approximately 1.5 Pa resolution). According to the manufacturer, the film used in this study is rated for normal pressures between 12 and 50 MPa, and has been validated independently by Selvadurai and Glaser (2015b) using a spherical indentation test. The pressure film is first placed between the fracture surfaces and compressed between the samples, in which it develops for a period of approximately 120 s. A specialized jig was constructed to ensure the repeatability of the pressure film measurements. Characterizing the fracture using the pressure film involved squeezing the film between the fracture surfaces under known uniaxial loads.

The pressure film offered a “one-time” estimate of the maximum normal stress — once a microcapsule becomes discolored, it does not return to its original state. For this reason, a new piece of pressure film is used to characterize the interface for each level of uniaxial loads. No shear stress (i.e., torque) was applied during compression because we are not interested in film-mediated sliding. The pressure film measures the same asperity configuration as our shear measurements, in which the equivalent uniaxial load is applied to press together the mated surfaces before any torque is applied.

After loading, the film was removed and optically digitized using an image scanner (MUSTEK ScanExpress A3 USB 2400 Pro Scanner). Algorithms in MATLAB were created to isolate, size, and catalog all contacting asperities in the static state. The light intensity of each microcapsule was converted to normal stress ($\sigma pix$) using the calibration curve documented in Selvadurai and Glaser (2015b). Once the image was converted to stress estimates, a lower threshold contact stress value was obtained in an iterative manner (Selvadurai and Glaser, 2012, 2015a). We assumed that the contact occurred only along the $r\u2010\theta $ plane (see Figure ^{11}) and that the force on each pixel was exactly perpendicular to the plane of the fracture (i.e., $z$-direction). For all pixels experiencing stresses above the threshold, a total reactive force $Fr$ was calculated (i.e., measured stress on pixel multiplied by the area of a pixel). The threshold was then varied iteratively until $Fr$ balanced the applied far-field uniaxial load, i.e., $Fr=Fn=\pi r02$$\sigma 0$. With the threshold determined, we were able to accurately measure the asperity contact area and normal stress, including mapping of the actual size, shape, and spatial distributions observed over the entire fracture surface (Selvadurai and Glaser, 2013).

### NUMERICAL MODEL

^{1}. The discontinuity in the body at $z=L/2$ is denoted as a positive (+) and negative (−) surface, which was necessary for defining the boundary conditions at this location. The entire model was composed of 76,200 eight-node linear brick elements (C3D8R), and it was refined more finely at the bonded/frictional interface (Figure

^{12}). Boundary conditions are as follows: