Estimation of the orientation of stress in the Earth’s crust without earthquake or borehole data

Mechanical stress acting in the Earth’s crust is a fundamental property that is important for a wide range of scientific and engineering applications. The orientation of maximum horizontal compressive stress can be estimated by inverting earthquake source mechanisms and measured directly from borehole-based measurements, but large regions of the continents have few or no observations. Here we present an approach to determine the orientation of maximum horizontal compressive stress by measuring stress-induced anisotropy of nonlinear susceptibility, which is the derivative of elastic modulus with respect to strain. Laboratory and Earth experiments show that nonlinear susceptibility is azimuthally dependent in an anisotropic stress field and is maximum in the orientation of maximum horizontal compressive stress. We observe this behavior in the Earth—in Oklahoma and New Mexico, U.S.A, where maximum nonlinear susceptibility coincides with the orientation of maximum horizontal compressive stress measured using traditional methods. Our measurements use empirical Green’s functions and solid-earth tides and can be applied at different temporal and spatial scales. The orientation of maximum horizontal compressive stress in the crust can be estimated in regions lacking borehole or earthquake source data through the use of ambient Earth noise to measure stress-induced anisotropy of nonlinear anelastic behavior.

K nowledge of the mechanical stress acting in the Earth's crust and lithosphere is important for a wide range of geophysical studies and applications [1][2][3][4] including plate tectonics 5,6 , seismicity and faulting [7][8][9][10][11] , and subsurface fluid behavior 9,12,13 . The stress field is commonly represented as the orientation of the maximum horizontal compressive stress (S Hmax ) 3,8,[14][15][16] , and other information regarding the principal components is often not known, much less the full stress tensor. At regional to tectonic-plate scales (100 s to 1000 s of km), the orientation of S Hmax is influenced by lateral plate boundary forces, tractions along the bottom of the lithosphere, and gravitational potential 17 . At local scales (<100 km), the orientation of S Hmax may vary owing to heterogeneities in density and elasticity, slipon faults 7,12 , and pore pressure 12 . The orientation of S Hmax is commonly estimated using borehole-based methods 18,19 , inverting earthquake focal mechanisms [20][21][22][23][24] , and less commonly by measuring the orientation of young, stress-sensitive geologic features 25 . Shear wave splitting and azimuthal seismic anisotropy are sometimes related to the stress field, but also reflect past deformation 26 . Borehole-based methods are high-cost, point measurements 27 , and are commonly applied in hydrocarbonproducing regions 8 . Interpreting earthquake focal mechanisms is limited to seismically active areas and requires an adequate monitoring network to produce high-quality, low-uncertainty source mechanisms. Because of limitations with these techniques, the stress field in broad regions of continental interiors is poorly constrained 15 .
Rocks are heterogeneous materials with stress and straindependent elastic properties, and finite, nonzero relaxation times (slow dynamics) [28][29][30] . This is in contrast to ideal linear elasticity (Hooke's Law), in which the elastic modulus is insensitive to strain, and the elastic response is instantaneous 31 . In individual rock samples, the relationship between stress, strain, and elasticity is complex 32,33 with mechanical damage and weak grain contacts being primarily responsible for nonlinear elastic behavior in brittle rocks 34 . Temperature, pressure, and the presence of fluids modulate the nonlinear behavior 34,35 . In the Earth, seismic velocities are commonly observed to be faster when rocks are compressed, usually interpreted as the closing of cracks and stiffening of internal contacts 36,37 , whereas they are typically slower after experiencing strong shaking, usually interpreted as the breaking or weakening of internal contacts 38,39 . After a dynamic disturbance, the material relaxes back to its original or a new metastable state via the process of slow dynamics 39,40 . Thus, rocks are metastable in their elastic behavior and strongly influenced by relatively weak external forces perturbing their material structure 35,38,39 .
Rock samples in laboratory experiments exhibit anisotropic linear and nonlinear elastic properties when differential stress is applied 41,42 . In two dimensions, differential stress is the difference between the two principal stresses. To demonstrate this nonlinear effect, Nur and Simmons (1969) applied uniaxial stress to a cylinder of granite, normal to the cylinder axis, and measured the travel time of an elastic wave as a function of angle with respect to the uniaxial stress 42 . The pressure derivative of the wave modulus with respect to strain (called nonlinear susceptibility, NS) was strongest when the angle between the uniaxial stress and the propagation of the probe wave was zero, and weakest when the angle was 90 degrees (Fig. 1). The effect was greatest for compressional P-waves, and the nonlinear elastic behavior can be quantified by measuring it. Stress-induced anisotropy in linear elasticity is exhibited by the velocity being faster for P-waves traveling in the same orientation of the applied uniaxial stress, rather than perpendicular. Stress-induced anisotropy in nonlinear elasticity is exhibited by the stress derivative of P-wave velocity being higher in the same orientation as the applied uniaxial stress, rather than perpendicular. Owing to a modest uniaxial or differential stress (1 MPa), one can observe anisotropy in linear elasticity of a few percent, whereas one can observe anisotropy in nonlinear elasticity of 100 s of percent 41,42 . Hence, it is preferable to measure anisotropy in nonlinear elasticity to infer principal stress orientations, rather than anisotropy in linear elasticity, as is done with shear wave splitting and azimuthal anisotropy of surface or body waves.
In addition to a nonlinear response to quasi-static stress, rocks exhibit creep behavior when compressed 43 . In this context, quasistatic means that the response of a system is as fast as the applied disturbance, and also that there is no inertia. Creep is delayed, time-dependent, and occurs when the strain response of a material is slower than the rate of the applied stress. The creep rate in rocks is highly sensitive to the magnitudes of the confining pressure and differential stress 43 . Yamamura et al., 29 measured elastic wave travel times in rock throughout the periodic tidal strain cycle and observed travel time differences from peak extension to the next peak compression (~6 h) were only weakly sensitive to the strain rate. Over several semi-diurnal (~12 h) cycles, the mean cycle travel time was most sensitive to the magnitude of peak extension, and relatively insensitive to the magnitude of peak compression. This indicates that the material response to extension was instantaneous while the material response to compression was delayed. The difference in travel time between peak extension and the next peak compression, and the mean cycle travel time, were strongly controlled by an apparently constant creep rate. The confining pressure and differential stress were constant in this study, but we expect the creep rate to be slower for higher confining stress and faster for higher differential stress, as it is in laboratory experiments on rocks 43 . Thus, under constant confining pressure, as is generally present in the Earth, we can use the creep rate to infer information about the differential stress.
We exploited this nonlinear elastic behavior in rocks and applied a technique to passively monitor the orientation of S Hmax in the lithosphere. Our approach to measure S Hmax orientation in situ relied on seismic velocity measurements that employed Empirical Green's Functions (EGF) derived from ambient Earth noise recorded at multiple pairs of seismic stations 44 . Crosscorrelations of diffuse seismic wavefields, such as ambient Earth noise or scattered coda waves, can be used to estimate Green's function. Most studies on Earth that measure temporal changes in seismic velocities do so by differencing the phase in the coda part of the EGF 36,38,39 . The coda of the EGF follows the direct waves and is the result of scattered waves that travel through some volume between the two stations 45 . We measured the velocity sensitivity to strain using a classic nonlinear acoustic approach known as the pump-probe method 33 , where, in a laboratory setting, the material is strained with a low-frequency oscillation (pump) and the elasticity is monitored by measuring the travel time of a high-frequency probe wave that was applied at different points in the pump cycle. For this study, solid-earth tides were used as the low-frequency pump and EGFs were the high-frequency probe. Solid-earth tides produce peak-to-peak axial strains up to 5 × 10 −8 46 , which corresponds to axial stress of~3.8 kPa using a P-wave velocity of 5 km s −1 for the shallow crust, corresponding to a reasonable bulk modulus of 63 GPa and Poisson's ratio of 0.25. The aerial strain tensor is elliptical with its long axis pointing towards the point on the Earth most directly facing the Moon (for Moon cycles) and Sun (for Sun cycles). For our study areas, there was an average ratio of~2.5 between the south-north and west-east axes of the strain tensor and this ratio is latitude-dependent. The magnitude of differential stress in Earth's brittle upper crust is uncertain in most areas, but is on the order of 10 s of MPa, decreasing toward the surface 11 . Consequently, the stress due to solid-earth tides is around three orders of magnitude less than the differential stress in the Earth's crust. By measuring changes in seismic velocity over the average tidal cycle, we could approximate the derivative of the seismic velocity with respect to strain, and also estimate the creep rate during compression. If the behavior we observe is similar to that shown by Yamamura et al. 29 , in a shallow cave, then the system in the Earth's upper crust is not quasi-static and creep behavior is important.
We performed this natural pump-probe experiment in two prototype studies located in north-central Oklahoma and northcentral New Mexico, USA (Fig. 2). We selected north-central Oklahoma because of the ongoing induced seismicity, generated by decades of injected wastewater from oil and gas operations 47,48 , that tends to occur on faults optimally orientated in the regional stress field 9 . North-central New Mexico was selected to test if we could resolve similar results in a geologic setting that straddles a continental rift separating the Colorado Plateau from the stable craton and has many different S Hmax orientations from Oklahoma 16,49 (Fig. 2). In north-central Oklahoma, S Hmax is oriented approximately N80E with some local variations 9 , but in north-central New Mexico, S Hmax is aligned nearly south-north along the Rio Grande Rift and rotates to a more east-west orientation eastward into northeastern New Mexico 16 . The dominant faulting style in Oklahoma is strike-slip, though there is some normal faulting in the northern part of our study area, whereas the faulting style is strongly normal faulting in northern New Mexico, associated with active extension along the Rio Grande Rift 8 .
Our results show that the Earth exhibits stress-induced anisotropy of NS that is aligned with S Hmax in these two different geologic settings, suggesting that the creep rate is fastest in the orientation of S Hmax . Since our measurements use only ambient seismic noise, there are several advantages over existing methods for estimating the orientation of S Hmax : (1) earthquake source properties are not used or required, (2) borehole measurements are not used or required, (3) sufficient seismic data exist in many regions of interest where traditional stress measurements are unavailable, and (4) the technique can be applied at a wide range of spatial and temporal scales.

Results
Our first measurement was to compare average seismic velocities in the Earth when tidal strains are extensional versus when tidal strains are compressional. We used the scattered wavefield in the coda of the EGFs calculated from the vertical components of the seismometers (Fig. 3). These scattered waves consisted of Rayleigh, compressional, and vertically polarized shear waves, though the relative amounts were difficult to discern. In both study areas, the Earth was slower during extension than during compression by fractional velocities of 0.04% with an uncertainty of ±0.003% for Oklahoma, and by 0.03% with an uncertainty of ±0.01% for New Mexico. This is consistent with softening of internal contacts during extension and stiffening of internal contacts during compression 29,36 . In Figs. 4 and 5, we report NS as fractional velocity change between velocities measured when the Earth is in compression to velocities measured when the Earth is in extension, though NS is actually fractional velocity change per unit strain. Tidal strain is on the order 10 −8 and it varies somewhat from cycle-to-cycle and by azimuth. Our results reflect the average peak-to-peak strain amplitude, which is discussed below.
Oklahoma. In Oklahoma, a sine function fitted to the results shows that the maximum (negative magnitude) NS occurs between 69 and 91°, depending upon the selection of stations (Figs. 2 and 4). Borehole measurements and a focal mechanism inversion show S Hmax orientations to be between 71°and 84°in the same region 9 . In that previous study, the reported S Hmax azimuth was 71°± 6°in the north (their area 2 N) and 82°± 6°in the south; (their areas 3 and 4, see Fig. 1 and Table 1 in Alt and Zoback 9 ). If we consider Fig. 4a-d versus Fig. 4e-f, we see a similar stress rotation from north to south. By considering S Hmax orientations from the 2016 release of the World Stress Map (WSM) 25 shown in Figs. 4g-l and 2, we can see that azimuths begin to decrease again (rotate more southwest-northeast) in the southernmost part of our study area. Some of the rotated stress indicators are outside our coverage area and we do not detect this rotation. So, this is a pointed disagreement between our measurements and the WSM. The areas considered in Alt and Zoback 9 , that we used for comparison, are not southward enough to include these more southwest-northeast oriented stress indicators. It is possible that we would detect this rotation with additional data to the south. Our measurements and these previously reported values reflect different spatial scales so exact comparisons are difficult. We note that considering only some of the northern stations produces ambiguous results because azimuthal variations are not clearly sinusoidal. This ambiguity and a few positive values in Fig. 4 that indicate velocities are faster during extension, maybe the result of poroelastic effects, and are discussed in more detail below. We estimated uncertainties in the predicted orientation of S Hmax by considering the uncertainties in the azimuthal measurements of NS (red bars on Fig. 4). The uncertainties in azimuth vary from 0.5 to 2.7 degrees depending upon the subarray. In addition, we used an f test to evaluate whether a sine function, which has three parameters (amplitude, phase, mean), is a statistically better model than a uniform function, which has one parameter (mean). Our results show that a sine function is a statistically better model for the observations with p-values ranging from 0.001 to 0.05, depending on the subarray. These results are reported in Fig. 4.
New Mexico. In New Mexico, a sine function fitted to the results from all nine stations (blue circles, Figs. 2 and 5) shows that the maximum (negative magnitude) NS occurs at 173° (Fig. 5). The maximum NS using only the six westernmost stations is 1°and using only the six easternmost stations is 161°, with the three central stations used in both subarrays. The regional published S Hmax orientations 9,25 show a transition, moving west to east, from slightly southwest-northeast to south-north within the Rio Grande Rift at this latitude (Fig. 2). Continuing east, a southeast-northwest S Hmax orientation may be expected if we interpolate between indicators in the Rio Grande Rift to indicators in northeastern New Mexico, although no published S Hmax orientations are available in the vicinity of the three eastern stations. Employing only the six western stations, the NS suggests  that S Hmax is 1°, which is consistent with S Hmax orientations observed within the rift valley and mountains to the west. Employing only the six easternmost stations, NS suggests that S Hmax is 161°, which is counterclockwise to the results of the six western stations, and intermediate between those reported within the Rio Grande Rift and the southeast-northwest S Hmax orientations reported east of our study area (east the blue circles in Fig. 2). Our NS-derived S Hmax results for all nine stations are in agreement with the average orientation of published S Hmax orientations within the footprint of the seismic array 9,25 . As no S Hmax orientations exist for the eastern part of the study area, our results using the eastern stations suggest the observed clockwise rotation from east to west transitions into our study area. The results provide a clear example of constraining the stress field using passive seismic data in a region where no other estimates are available. We estimated uncertainties in the predicted orientation of S Hmax , and tested a sine function versus a uniform function as models to describe the azimuthal variations, in the same manner as the Oklahoma results. The uncertainties in azimuth vary from 3.0 to 4.4 degrees depending upon the subarray. Our results show that a sine function is a statistically better model for the observations when using all nine stations with a p value of 0.0005. Using only the six western stations, the p value is 0.05, and using only the six eastern stations the p value is 0.5. So, for the eastern stations, we cannot reject the null hypothesis that a uniform function describes the observation better than a sine function. These results are reported in Fig. 5.
The uncertainties that we report on Figs. 4 and 5 represent 1-sigma standard deviation measurement uncertainties for average stress orientations and not the variability in stress orientations within the corresponding seismic array. It is possible to have low uncertainty for average behavior from a set of individual measurements with higher uncertainties. To explore this further, we can consider four panels, Fig. 4a-d. As we remove northern stations, the remaining interstation paths produce a lower azimuth (74°-69°), then higher (69°-74°), then lower again (74°-69°) with measurement uncertainties of~1°. This suggests something like a striped pattern in the spatial distribution of stress orientations. This could be true, but it also leads us to believe that our uncertainties may be underestimated by up to a few degrees, and perhaps the true average for this area is just  somewhere between 69°and 74°. The contribution of data and measurement uncertainties from each pair of stations, to the final result, is not linear as with simple averaging.

Discussion
Here, we discuss different mechanisms that may explain our observation that the fractional velocity change (Δv v −1 ) between tidal extension and tidal compression (Δv Δε −1 ) was the greatest in the orientation of S Hmax . In general, we assumed that measured velocity variations are related to the stiffness across internal contacts 41,42 . These contacts, such as grain boundaries or fractures, produce nonlinear behavior and have a direction in which they open (weaken) and close (strengthen) 28,50,51 . We assumed that internal contacts naturally exist at all orientations in rocks of the subsurface, and considered only those that open and close in a horizontal direction because these internal contacts are most strongly affected by the differential, areal stress. From this population, we considered how internal contacts with different orientations respond to tidal forcing and the ambient stress field (S Hmax ) together. When we refer to contact orientations below, we mean the orientation in which they open and close. For both P-waves and S-waves, velocities are sensitive to contacts oriented in the same direction as wave propagation, but whose sensitivities are higher for P-waves than for S-waves 41,42 . Both Love and Rayleigh waves would be considered SH waves as described in Nur and Simmons 42 and Johnson and Rasolofosaon 41 because they consider a vertically applied uniaxial stress, and in the Earth the uniaxial (differential) stress is horizontal. A shear wave would have to travel vertically and be polarized in the orientation of the uniaxial stress to be considered SV in their notation. Because we use vertical component data, we expect that our wavefield consists of scattered Rayleigh waves, Pwaves, and V sv (vertically polarized) shear waves. According to laboratory and theoretical results, the pressure derivative of Pwaves, at small pressures, is~4× greater for waves traveling in the same direction as the applied uniaxial stress than for waves traveling perpendicular to the applied uniaxial stress, and~2× greater for SH waves 41,42 . Even though SH waves are much less sensitive to differential stress than P-waves in this regard, these amounts are still dramatically higher than expected anisotropy in the wave speeds themselves, which would be only a few percent for differential stress of a few MPa 41,42 .
To determine our preferred mechanism, presented first, we used the results from Yamamura et al. 29 , and laboratory experiments as guidance. In the Yamamura study, Δv v −1 over one tidal cycle was only weakly sensitive to the strain rate, and there was no unique correspondence between velocity and strain, suggesting that Δv v −1 exhibits creep during compression 29 . We knew from laboratory experiments on rocks that creep rate is highly sensitive to confining pressure and differential stress 43 . Rocks at any given depth in the Earth have the same overburden with no azimuthal dependence, so overburden could not produce an anisotropic result. Spatial variations in pore pressure at the same depth could result in different effective confining pressures, which would affect creep rates, but also would not produce an anisotropic result unless poroelastic properties were highly anisotropic.
The data from Yamamura et al. 29 were measured at a single azimuth and the orientation of S Hmax was unknown to us, so no azimuthal dependence could be discerned. However, we considered the effect of the differential ambient stress field, which is three orders of magnitude stronger than the tidal stress (~10 MPa versus~4 kPa). For internal contacts in the orientation of S Hmax , the differential stress always promotes strengthening or closure and results in faster creep rates for closing contacts during compressive tidal stress 43 . For internal contacts perpendicular to S Hmax (i.e., in the orientation of S Hmin , the minimum horizontal compressive stress), the differential stress always opposes closure and results in slower creep rates. These contacts experience negative (contact weakening) differential stress even when tidal stress is compressive because the magnitude of the ambient stress field is three orders of magnitude greater. Therefore, the creep rate during tidal compression should always be faster for contacts in the orientation of S Hmax than for contacts in the orientation of S Hmin . Since in Yamamura et al. 29 , Δv v −1 was determined primarily by the creep rate, we should expect to observe the highest Δv v −1 values in the orientation of S Hmax .
The Yamamura et al. 29 experiment was at very low confining stress in soft rock; P-wave velocity~2 km s −1 , saturated, and with 40% porosity. The rocks in our study were likely less porous, at least 2-3 times stiffer, and were likely to be saturated 52 . Both the confining pressure and the differential stress were several orders of magnitude higher. However, creep behavior is ubiquitous in rocks in the laboratory over a wide variety of types and conditions 43 .
A second possible mechanism is that ambient tectonic stress and tidal stress are combined in a nonlinear manner. In the laboratory experiment described previously, an increasing uniaxial stress-induced elastic anisotropy in dry granite without any confining pressure 42 . These results were for a quasi-static system and did not consider behavior during loading and unloading separately 42 . Our experiment consisted of constant confining and differential stresses under likely wet conditions, combined with a near-isotropic cyclical (tidal) stress. A simple application of the laboratory results would be that the NS decreases with increasing magnitude of the uniaxial stress in any specific orientation, which would produce the opposite result to the one we obtained. However, one should not assume quasi-static conditions in our study. A more applicable experiment to our study would provide a better comparison.
A third possible mechanism for the observed azimuthal dependence of fractional velocity change in the study areas is that there is a preferred orientation for cracks, internal contacts, or anisotropic minerals that are not related to the present-day stress field 53 . Since we were not measuring anisotropy in linear elasticity (Δv v −1 ), but instead anisotropy in nonlinear elasticity (Δv Δε −1 ), it is not clear that the same behaviors observed in linear elasticity, such as shear wave splitting or azimuthal anisotropy in surface waves, apply to nonlinear elasticity. If a material were merely softer in one orientation than another, its modulus is not necessarily more or less sensitive to strain.
A fourth possible mechanism is that apparent azimuthal variations in the measured NS (Δv Δε −1 ) were the result of biases in the ambient noise field 54 . An intrinsic assumption when using cross-correlations to recover EGFs was that ambient noise was an equipartitioned, diffuse wavefield 55 . Any non-equipartitioned, non-diffuse properties of the ambient wavefield would introduce biases that could give incorrect velocity measurements. Though undoubtedly, there were non-equipartitioned, non-diffuse aspects to the ambient wavefield, even after pre-processing, we do not think it was sufficient to invalidate our results. The data from the two study sites were not recorded at the same time, but they would have similar ocean-based sources and seasonal variations 55 and yet, they provided different, and apparently correct estimates of S Hmax . Also, as we were measuring velocity changes, owing to solid-earth tides, rather than absolute velocities, consistent or slowly-changing biases would not affect our results.
The difference in strain between maximum extension and maximum compression for solid-earth tides was of the order 5 × 10 −8 according to model calculations made with the software package SPOTL 46 . As we observed fractional changes in the velocity of~0.04%, associated with these differences in tidal strain, the observed NS was of the order 10 4 . Our results are similar in magnitude to those found by Takano et al. 56 near a volcano in Japan using EGF frequencies of 1-2 Hz, and an order of magnitude higher than those found by Hillers et al. 57 in California using frequencies of 2-8 Hz. It is difficult to compare these results to ours beyond the raw numbers because the apertures of the arrays in the previous studies were both under 1 km, and ours were >100 km. Owing to the larger apertures and lower frequencies in our study, we are likely probing deeper into the subsurface and are certainly averaging over a larger area. The previous studies measured nonlinearity, but did not report azimuthal differences, nor the relationship between nonlinearity and stress-induced anisotropy.
In addition to the softening and stiffening of internal contacts, there may be poroelastic effects owing to interactions between solid-earth tides and pore pressure. In saturated conditions, pore pressure increases during applied compression and decreases during applied extension, with pore pressure having the opposite effect on the effective confining stress as the applied tidal stress. If all station pairs in an array experience the same poroelastic conditions, the effect of stress-induced anisotropy is preserved, though curves shown in Figs. 4 and 5 would shift upward (positive). In some cases, we might even observe positive values 57 . If different station pairs experience different poroelastic conditions, such as different coupling between pore pressure and effective stress, then the effect of stress-induced anisotropy may not be as evident. In Oklahoma, we were able to get good estimates for S Hmax despite possible contributions from heterogeneous poroelastic conditions 58 that were apparent when using fewer station pairs. This may be why using only northern stations in Oklahoma did not produce a sinusoidal pattern with azimuth, and that the sinusoidal fit was generally worse when using fewer station pairs in Figs. 4 and 5.
Measuring and modeling principal stress orientations in the Earth's crust is challenging, and it is important to match the length scale and depth to the desired application. Since stress heterogeneity likely exists at all scales 10,20,27 , it is advantageous to measure S Hmax at the same length scales and depths as the application. EGFs can be calculated using varying interstation distances to provide S Hmax estimates at different horizontal length scales. Estimating S Hmax at specific depths is more challenging but possible when using relative depths inferred through frequency content and coda time offset in the EGFs 38,45 . Perhaps most importantly, measuring the orientation of S Hmax is not limited to locations with earthquakes or boreholes, and provides data-driven constraints to regional estimates. Additional possibilities include calculating the time evolution of NS, which could reveal temporal changes in relative amplitude and orientation of S Hmax . Temporal changes may occur in the immediate vicinity of a recent earthquake, on a plate boundary experiencing stress/strain loading, or in a reservoir being depleted.
Stress and deformation patterns in tectonic plates are generated from global and local scale mantle convection 59,60 , gravitational potential 61 , and plate boundary tractions 16,62 . The relative contribution of these mechanisms is unknown. Ultimately, we do not know to what extent continental scale stress models represent the actual stress field in regions with few or no measurements to constrain these estimates. Therefore, we cannot attempt to model or characterize mechanisms for these unknown heterogeneities. This method makes possible dense and uniform observations of the orientation of S Hmax across continental regions, which will improve stress models and thereby our understanding of the underlying geodynamical processes. One possible application would be to apply this method at a continental scale by using many overlapping subarrays within Earthscope's US Array, similar to what we have done in New Mexico. Another option would be to design a denser, local array over a particular area of interest, such as a geothermal field.
In this study, we calculated EGFs as a function of tidal strain and azimuth in north-central Oklahoma and north-central New Mexico to constrain NS and estimated the orientation of S Hmax . Our results in both study areas show that seismic velocities were, on average, faster when tidal strains were in compression relative to when they were an extension. We observed stress-induced anisotropy in nonlinear anelastic behavior, which is aligned with S Hmax and provides a technique to estimate the orientation of S Hmax without focal mechanism inversions or borehole measurements. Large-scale application of this method may resolve additional tensor properties of the nonlinear behavior, reveal how S Hmax varies with horizontal length scales and depth, and how S Hmax evolves temporally in areas such as fluid reservoirs and active fault zones.
Methods and data Data. We used publicly available seismic data from the two study areas, north-central Oklahoma and north-central New Mexico. For Oklahoma, we obtained waveform data recorded by the Nanometrics Research Array (NX) from the Incorporated Research Institutions for Seismology Data Management Center 63 . The NX array consists of 30 broadband, three-component instruments that recorded 100 samples per second for~3 years between mid-2013 and mid-2016 (red circles in Fig. 2). For New Mexico, we obtained waveform data from nine stations in Earthscope's Transportable Array (TA) also from the Incorporated Research Institutions for Seismology Data Management Center 64 . This subarray consists of nine broadband, threecomponent instruments that recorded 40 samples per second for~2 years between mid-2008 and mid-2010 (blue circles in Fig. 2). We used only the vertical component for all seismic data.
Signal processing. Signal processing steps were performed using the Python package ObsPy 65 . We organized the data into daylong segments and deconvolved the instrument response. When calculating EGFs from continuous broadband seismic data it was important to remove transient signals like earthquakes 66 . We removed earthquake signals from the data using earthquakes identified in the U.S. Geological Survey Comprehensive Catalog. We multiplied windows within the waveforms by zero, tapering to 1 at the beginning and end, following three sets of earthquake criteria: (1) earthquakes with a minimum magnitude of 3.5 and a maximum distance of 30 km from the array, between surface wave velocities of 2 and 5 km s −1 , (2) earthquakes with a minimum magnitude 5 and a maximum distance of 2000 km from the array, between surface wave velocities of 2 and 7 km s −1 , and (3) earthquakes with a minimum magnitude of 6 at any distance, between surface wave velocities of 2 and 8 km s −1 . These windows were chosen because they bracket the expected start and end of the wave train from the corresponding earthquakes. There may be smaller earthquakes present on the recordings, but their shorter durations and weaker amplitudes are less of a concern for interfering with our measurements. This resulted in the zeroing of 7.9% of waveforms for Oklahoma and 4.3% of waveforms for New Mexico. The disparity exists because there were more local earthquakes in Oklahoma than in New Mexico during the study intervals. Additionally, we clipped all signals greater than three times the RMS of each day-long segment to remove nonearthquake signals observed as emergent or impulsive noise.
Tidal strain. In order to determine time windows that were in periods of high (extensional) and low (compressional) strain, we determined the volumetric tidal strain using the software package SPOTL 46 . We used the volumetric strain component because we expect the nonlinear behavior to be localized on pre-existing internal contacts, which may have any orientation. We divided time into two groups according to tidal strain magnitude, the top 25% and the bottom 25%, where top refers to maximum extension and the bottom refers to maximum compression. For stress modeling in the Earth, a common convention is that compressive stress has a positive sign, because absolute stress is nearly always compressive. Except when referring to S Hmax (maximum horizontal compressive stress), we use the opposite convention so that positive axial stress results in extension if the elastic modulus is also positive. In this manuscript, we frequently use the terminology compressive and extensional to avoid any misunderstandings.
Empirical Green's functions. We cut and merged the preprocessed waveforms for all stations in order to separate the data into two groups corresponding to time periods when tidal strains were in compression and time periods when tidal strains were in extension. We discarded any waveform segments shorter than 30 min duration because very short segments do not contain enough data to be useful. We calculated an EGF for each selected station pair, for both compressional and extensional groups using a phase cross-correlation method 67 in which we pre-whitened the spectrum before applying a phase cross-correlation. There were~780 individual segments for each station pair during the recording period for both study areas, though the actual number for each pair varied based on data availability, data quality, and other factors. In particular, we discarded any segments that produced a velocity measurement that exceeded four times the standard deviation from the mean of all the segments. For the Oklahoma stations, we empirically determined that a station separation distance between 30 and 60 km produced the best EGFs and selected station pairs accordingly. This involved applying a bandpass filter with corners of 0.1 and 1 Hz, plotting the EGFs, and visually inspecting them. In our visual inspection, we looked for a well-defined, dispersive Rayleigh wave with decaying amplitude in the coda, and when plotting the EGFs by interstation distance, there was a clear and consistent move-out of the waveforms. For the New Mexico stations, we used all pairs because there were fewer stations.
Next, we describe the EGF stacking procedure (Fig. 3). We selected 14-day windows, and for each station pair, we selected all EGFs whose segment start time fell within ±7 days of the center of the window. Each EGF was scaled by the square root of the duration of the underlying time series before stacking so that each EGF contributed to the stack according to the amount of data it contained. We calculated the Pearson correlation coefficient 68 , which varies between −1 (perfectly anti-correlated) and 1 (perfectly correlated), of each EGF with the stack. Any EGF yielding a Pearson correlation coefficient <0.5 was discarded and then we produced a new stack with the remaining EGFs. We reevaluated the discarded EGFs and any that had a Pearson correlation coefficient >0.5 using the updated stack was reinserted to create a new stack. The process was repeated until there were no discarded EGFs with a Pearson correlation coefficient greater than 0.5 with the stack. Subsequent 14-day windows were calculated with a 7-day overlap. This stacking procedure was intended to include as many observations as possible while discarding outliers. The process resulted in stacked EGFs for each station pair that represented 14-day windows with 7-day overlaps for each of the two strain bins described above (Fig. 3). We estimated the measurement uncertainty using the standard deviation of the 14-day windows.
We summed the positive and negative lagging parts of the EGF and selected the coda part as shown in Fig. 3 to avoid direct wave arrivals. We determined the average phase difference and velocity change (Δv v −1 ) between the two EGFs (compression and extension) in a 30 s coda window for waves between 4 and 5 s period following the steps outlined in the wavelet method of Mao et al. 69 . This method uses a continuous wavelet transform to convert time series data (the coda) into the frequency domain. To perform the continuous wavelet transforms, we used the Python package PyWavelets 70 . Unlike a Fourier transform, a wavelet transform is localized in both time and frequency because the wavelet has both a finite time-width and frequency bandwidth. We selected periods of 4-5 s because Rayleigh waves in this period range are sensitive to structures in the top few km of the Earth. The coda also contained scattered body waves consisting of P and vertically polarized S waves (V sv ).
Once the coda waveforms were represented in the frequency domain, we could measure small phase shifts as a function of frequency and time offset in the coda. We selected a coda window of 30 s because that is the period of time when we expect scattered wave arrivals based on past experience 38 . We used a Morlet wavelet 71 with ω 0 = 0.25 Hz that corresponded to the periods we analyzed and allowed us to recover the known phase shifts in simple synthetic examples. In these synthetic examples, we took a real EGF coda waveform and manually introduced phase shifts of various magnitudes and signs by stretching or compressing the waveform in time. Then, we verified that we could recover the apparent velocity change associated with the imposed phase shift. This method had the ability to measure phase shifts associated with changes in velocity as a function of frequency and coda offset time. In general, the earlier part of the coda contained more scattered surface waves, whereas the latter part contained more body waves, with the transition time governed by the scattering properties of the subsurface 45 . Scattered Rayleigh waves are sensitive to the upper 2-3 km for periods of 4-5 s. If the window of the measurements contained body waves (P and V sv ), then the waves would be sensitive over a greater depth range, depending on the velocity of the scattered waves, and the scattering properties in the subsurface. We were only interested in the average velocity changes for this analysis, and so we calculated the average phase shift, and associated velocity change, over the entire coda window and period range. Results measured as a function of frequency and coda offset time likely contain depth information 38,45,69 but were beyond the scope of this study. Along with measuring phase shifts, we also calculated coherence between EGFs and discarded any cases where the average coherence fell below 0.95. In our convention, a negative (Δv v −1 ) value meant that according to the EGF coda on the vertical channel, the Earth was slower during extension than during compression.
Azimuthal measurements. We grouped and stacked the station pairs by azimuth so that we could examine any orientation dependence on the results. The azimuth for each pair was determined using the relationship from the more western station to the more eastern station so that values are always between 0 and 180 degrees, with 0 and 180 indicating south-north and 90 indicating west-east. For Oklahoma, we considered nine azimuths in 20-degree steps. For each azimuthal interval, we averaged the Δv v −1 values for all pairs whose azimuth was within ±20 degrees with wrapping. For example, at an azimuth of 0 degree, we averaged paths between 0 and 20 plus those between 160 and 180 degrees (which is equivalent to between −20 and 0 degree). For New Mexico, we considered all station pairs individually because there were not enough pairs to average in azimuthal bins. In addition to calculating the average Δv v −1 at different azimuths, we fit a sine function, periodic on 2θ, to the results. To obtain the predicted values for S Hmax reported on Figs. 4a-f and 5a-c, we generated 1000 realizations where we drew a value for fractional velocity change for each station pair (New Mexico) or azimuthal bin (Oklahoma) from a normal distribution using the mean and standard deviation indicated by the red bars, then solved for the best fitting sine function. We report the mean and standard deviation of the 1000 realizations on Figs. 4a-f and 5a-c. Generating more realizations does not change the result. We also applied an f test to evaluate whether or not a sine function is a statistically better model than a uniform function to describe the azimuthal observations. The sine function has three parameters (amplitude, phase, mean), whereas the uniform function has only one parameter (mean). Applying the appropriate penalty for having more model parameters, a low p value indicates that the sine function is a better model to describe the observations. The p value of the f test is annotated in Figs. 4a-f and 5a-c. Only the results for the eastern stations in New Mexico indicate that a sine function is not a statistically significantly better model than a uniform function.
Uncertainties in our velocity measurements were difficult to precisely estimate. When calculating EGFs, there was an intrinsic assumption that noise sources are equipartitioned, white, azimuthally uniform, and stationary. This is never true on Earth, but we could take steps to reduce the influence of recordings that violate these assumptions. We assumed that EGF coda velocities in our study area vary by no more than the amounts observed in other studies, <±1% 39,56,57 . As we never directly compared data that were recorded >14 days apart, seasonal variations were not expected to be important, including spectral content, azimuthal variations, and relative amounts of coherent and incoherent noise. By discarding EGFs that were not well correlated and coda that was not highly coherent, we avoided time segments that were badly contaminated with non-equipartitioned waveforms. In both study areas, we stacked over two years of data. The resulting uncertainty based on the variance in the accepted measurements suggested that uncertainties are no more than one-tenth (Oklahoma) and one-third (New Mexico) of the magnitude of the measurements for the full dataset. Uncertainties for individual stations (New Mexico) and azimuthal bins (Oklahoma) are indicated in Figs. 4 and 5. Average uncertainties on the full dataset for the two study areas include the fact that the Δv v −1 magnitudes are different at different azimuths. Therefore, these values overestimate the actual measurement error in this regard.

Data availability
All data used in this study are available through the Incorporated Institutions for Seismology Data Management Center (www.iris.edu). Information regarding the Nanometrics Research Array can be found here (ds.iris.edu/mda/NX/) and Earthscope's Transportable Array here (ds.iris.edu/mda/TA/). Time series data can be requested through several tools provided here (ds.iris.edu/ds/nodes/dmc/tools/#data_types = time series).

Code availability
This study was performed using the Python package Obspy 65 and PyWavelets 70 and uses workflows provided in Ventosa et al. 67 and Mao et al. 69 .