Measuring S Hmax with Stress-Induced Anisotropy in Nonlinear Anelastic Behavior

Mechanical stress acting in the Earth`s crust is a fundamental property that has a wide 8 range of geophysical applications, from tectonic movements to energy production. The 9 orientation of maximum horizontal compressive stress, S Hmax can be estimated by inverting 10 earthquake source mechanisms and directly from borehole-based measurements, but large 11 regions of the continents have few or no observations. Available observations often represent a 12 variety of length scales and depths, and can be difficult to reconcile. Here we present a new 13 approach to determine S Hmax by measuring stress induced anisotropy of nonlinear susceptibility. We observe that nonlinear susceptibility is azimuthally dependent in the Earth and maximum 15 when parallel to S Hmax , as predicted by laboratory experiments. Our measurements use empirical 16 Green’s functions that are applicable for different temporal and spatial scales. The method can 17 quantify the orientation of S Hmax in regions where no measurements exist today.


Introduction 20
Knowledge of the mechanical stress acting in the Earth's crust and lithosphere is important 21 for a wide range of geophysical studies and applications [1][2][3][4] including plate tectonics 5,6 , seismicity 7-22 11 and subsurface fluid behavior 9,12,13 . It is commonly represented as the orientation of the 23 maximum horizontal compressive stress (SHmax) 3,8,14-16 . Other information regarding the principle 24 components is often not known, much less the full stress tensor. At regional to tectonic-plate 25 scales, the orientation of SHmax is determined by plate boundary forces and tractions along the 26 bottom of the lithosphere 17 . At local scales (<10 km), the orientation of SHmax may vary due to 27 heterogeneities in density and elasticity, slip on faults 7,12 , and pore pressure 12 . The orientation of 28 SHmax is commonly estimated using borehole-based methods 18,19 and inverting earthquake focal 29 mechanisms 20-24 , and less commonly by measuring the orientation of stress sensitive geologic 30 features 25 .
Borehole-based methods are high-cost, point measurements with unknown 31 applicability away from the borehole 26 , and commonly applied in hydrocarbon producing regions 8 . 32 Interpreting earthquake focal mechanisms is limited to seismically-active areas and requires an 33 adequate monitoring network. Because of the limitations of these techniques, broad regions of 34 continental interiors are poorly constrained where there are few or no measurements 15 . 35 Rock samples in laboratory experiments typically exhibit anisotropic nonlinear elastic 36 properties when a uniaxial stress is applied 27,28 (Figure 1). In the laboratory, the pressure derivative 37 of the wave modulus (called nonlinear susceptibility, NS) is strongest when the angle between the 38 uniaxial stress and the propagation of the probe wave is zero, and weakest when the angle is 90 39 degrees. The effect is greatest for compressional P-waves and the nonlinear elastic behavior is 40 stations, the NS predicts SHmax at 8°, which is consistent with stress indicators within the rift 115 valley and mountains to the west. Using the 6 easternmost stations, NS predicts SHmax at 163°, 116 which is rotated counterclockwise from the results of the western stations, and intermediate 117 between the reported stress indicators within the Rio Grande Rift and the reported southeast-118 northwest stress indicators east of the study area. Our NS derived SHmax results for all 9 stations 119 are in agreement with the average orientation of stress indicates within the footprint of the 120 seismic array. Since no stress indicators exist for the eastern part of the study area, our results 121 using the eastern stations suggest the observed clockwise rotation from east to west transitions 122 into our study area. The results provide a clear example of constraining the stress field using 123 passive seismic data in a region where no other estimates are available. 124

Discussion 125
The azimuthal dependence of NS closely tracks the orientation of SHmax in the Earth as 126 shown in the laboratory 27,28 . In terms of elastic constants, what we measure is closely related to 127 the 1-D nonlinear anelastic coefficient β which is the coefficient linearly related to strain in a 128 Taylor expansion of Hooke's law (see, e.g., equations 7 and 8 in Johnson and Rasolofosaon 27 129 and Pantea et al. 45 ). In single crystals and metals β is less than 10. In Earth materials it can be 130 considerably larger, order 10 2 -10 3 underscoring how very nonlinear elastic Earth materials can 131 be. In the laboratory experiments, stress induced anisotropy is strongest for P-waves 28 , which 132 suggests our measurements describe the scattered compressional-wave energy either in the form 133 of Rayleigh waves or P-body waves. Since our measurements are made with the vertical 134 channels of a seismic station pair, we are measuring a specific component of the nonlinear 135 coefficient β, which may be better represented as a tensor. Exploring the possible tensor 136 properties of β is beyond the scope of this study but may be possible by developing a 6 component NS tensor using all combination of station channels. Such a measurement would be 138 very useful in discerning the anisotropic mechanical damage variations in the upper crust. 139 The difference in strain between maximum extension and maximum compression for 140 solid Earth tides is of the order 5x10 -8 , which means the observed NS (dv/dε, change in velocity 141 over change in strain) is of the order 10 4 -10 5 . These results are similar in magnitude to those 142 found by Takano et al. 46 near a volcano in Japan using EGF frequencies of 1-2 Hz, and an order 143 of magnitude higher than those found by Hillers et al. 47 in California using frequencies of 2-8 144 Hz. In these two previous studies, the authors measured nonlinearity, but did not report 145 azimuthal differences, nor the relationship between nonlinearity and stress-induced anisotropy. Measuring and modeling stress in the Earth's crust is challenging and it is important to 161 match the length scale and depth to the desired application. Since stress heterogeneity likely 162 exists at all scales 10,20,26 , a measured stress or SHmax may not be representative of different length 163 scales or depths. Our method has the potential to address some of these challenges. EGFs can 164 be calculated using varying interstation distances to provide SHmax estimates at different 165 horizontal length scales. Estimating SHmax at specific depths is more challenging but possible 166 when using relative depths inferred through frequency content and coda time offset in the 167 non-earthquake signals observed as emergent or impulsive noise. 217

Tidal Strain 218
The volumetric tidal strain is obtained using the software package SPOTL 54 . We use the 219 volumetric strain component because we expect the nonlinear behavior to be localized on pre-220 existing faults, which may be at any orientation. We divided time into segments that fit in two 221 strain magnitude bins, the top 25% and the bottom 25% where "top" refers to maximum extension 222 and "bottom" refers to maximum compression. 223

Empirical Green's Functions 224
We cut and merge the day long preprocessed waveforms into segments appropriate for 225 each stress bin and discard anything shorter than 30 minutes. We empirically determined that a 226 station separation distance between 30 and 60 km produce the best EGFs for Oklahoma and 227 selected station pairs accordingly. For the New Mexico stations we used all pairs. We 228 calculated an EGF for each selected station pair, and segment for all bins using a phase cross 229 correlation method 55 where we pre-whiten the spectrum before applying a phase cross 230 correlation. There are approximately 780 segments for each station pair during the recording 231 period, though the actual number for each pair varies based on data availability, data quality, and 232 other factors. 233 Next, we describe the EGF stacking procedure ( Figure 5). For each station pair we selected 234 14 day windows, and using the center of the window, selected all EGFs whose segment start time 235 falls within +/-7 days. Each EGF is scaled by the square root of the duration of the underlying 236 time series and are all stacked. We calculate the Pearson correlation coefficient of each EGF with 237 the stack. Any EGF that has a value of less than 0.5 is discarded and we produce a new stack with 238 the remaining EGFs. We reevaluate the discarded EGFs and any that have a Pearson correlation coefficient greater 0.5 using the updated stack is re-included to create a new stack. The process is 240 repeated until there are no discarded EGFs with a Pearson correlation coefficient greater than 0.5. 241 Subsequent 14 day windows are calculated with a 7 day overlap. This stacking procedure is 242 intended to include as many observations as possible while discarding outliers. The outcome is 243 stacked EGFs for each station pair that represent 14 day windows with 7 day overlap for each of 244 the two strain bins described above. 245 We sum the causal and acausal parts of the EGF and select the coda part as shown in Figure  246 5 to avoid direct wave arrivals. We determine the average phase difference and velocity change 247 Along with measuring phase shifts, we also calculate coherence between the two EGF stacks and 255 discard any cases where the average coherence falls below 0.95. In our convention a negative 256 (Δv/v) value means that the Earth is slower during extension than during compression, which we 257 EGFs, there is an intrinsic assumption that noise sources are equipartitioned, white, azimuthally 278 uniform, and stationary. This is never true in the Earth, but we can take steps to reduce the 279 influence of recordings that violate these assumptions. We assume that velocities in our study 280 area vary measurably by no more than the amounts observed in other studies, less than +/-281 1% 39,46,47 . Since we never directly compare data that has been recorded more than 14 days apart, 282 seasonal variations are not expected to be important, including spectral content, azimuthal 283 variations, and relative amounts of coherent and incoherent noise. By discarding EGFs that are 284 not well-correlated and coda that are not highly coherent we avoid noise or signals that are not 285 stationary. In both study areas, we stack over two years of data. The resulting uncertainty based on the variance in the accepted measurements suggests that uncertainties are at least an order of 287 magnitude less than the magnitude of the measurements. 288