Shear-Wave Anisotropy Measurements in the Crust from Receiver Functions: An Interplay of Lower and Upper Crustal Anisotropy

Abstract

We report a study using teleseismic P-wave receiver functions to infer the orientation of the maximum horizontal principal stress from the direction of upper crustal shear-wave velocity anisotropy. We apply an inverse approach using the Neighborhood Algorithm to conduct a nonlinear search, attaining a best-fitting crustal model that includes shear velocity anisotropy. Unlike previous methods reported in the literature, this method is able to distinguish anisotropy in the upper, brittle crust from that in the lower, ductile crust in certain instances. We apply this method to teleseismically recorded earthquakes in the Central Valley of California, the Permian Basin, Texas, northern Oklahoma and sites near the San Andreas Fault in California. Of the forty-one stations to which we apply this method, twenty have a good apparent signal. A misfit calculation is performed by calculating a zero-lag cross-correlation coefficient for each modeled receiver function with the data for a given back azimuth range. While the fast polarization direction in the upper crust of some of these stations aligns with independent indicators of the direction of the maximum horizontal principal stress, the fast direction in the upper crust at other stations does not, apparently indicating that the anisotropy was resulting from a different mechanism.

1. Introduction

Knowledge of stress orientation and relative magnitude in the crust are important for geo-energy development and geologic carbon storage. These parameters generally come from four sources of information: wellbore data (principally observations of stress-induced compressional and tensile failures in relatively deep wellbores), earthquake focal plane mechanisms (ideally the inversion of multiple high-quality mechanisms in a given area), relatively recent geologic indicators (volcanic alignments or observations of fault slip) and shear velocity anisotropy from local earthquakes or anthropogenic sources [1]. As originally demonstrated by [2], the application of stringent quality criteria results in remarkable consistency among these methods despite the fact that the depths and volumes of rock sampled vary markedly. Over the past 40 years, there has been considerable improvement in stress determination methodologies as well as the amount of stress data available. Nonetheless, considerable gaps in coverage persist both at local and regional scales in North America and other areas around the world [3,4,5]. In this paper, we present a study investigating whether crustal shear velocity anisotropy determined from receiver functions of teleseismic earthquakes could provide information about stress orientation where no other stress information is available.

Shear-wave velocity anisotropy in the crust could result from two processes—the anisotropy could be related to the current stress state or there could be an intrinsic structural anisotropy due to aligned fractures, fabrics, or foliation. In the former case where the tectonic stresses are controlling the anisotropy, if the upper, brittle portion of the crust contains pre-existing fractures and faults at many orientations, the faults and fractures oriented roughly perpendicular to the direction of the maximum horizontal principal stress (S_Hmax) will close and/or there could be dilation of faults and fractures oriented parallel to S_Hmax. Both processes cause the fast polarization direction of the shear-wave anisotropy to align to the direction of S_Hmax [6,7]. The degree of anisotropy associated with the closure and dilation of faults and fractures is on the order of 2–4% [8,9]. The latter phenomenon that causes shear-wave anisotropy is the structural alignment of rock fabrics, which is typically greater in the lower crust where ductile flow can align minerals [10]. The degree of anisotropy associated with structural alignment can be as high as ~20% [11]. Due to the difference in the magnitudes of anisotropy from the two phenomena, we expect that if the contribution from fabric-induced anisotropy in the upper crust is slight, then it should be possible to determine the direction of S_Hmax from the fast polarization direction. On the other hand, if there is significant structurally controlled anisotropy, especially in the lower crust where the effect is amplified by the ductility, we show that the ability to determine the direction of S_Hmax can be complicated.

We test the inversion to perform a nonlinear search of parameter space to minimize the misfit between synthetic data and the recorded receiver function data of forty-one stations. The remainder of this paper details the methodology of this approach and then analyzes the results to examine the interplay between lower and upper crustal anisotropy and the potential for measuring S_Hmax in such a manner.

2. Materials and Methods

Forty-one stations were analyzed using the Neighborhood Algorithm inversion [12,13]: four in the Central Valley of California, four in southern California, seventeen in the Permian Basin, fifteen in Oklahoma and Kansas, and also the Parkfield station overlying the San Andreas Fault in central California (Figure 1). The southern California stations were chosen because local earthquake shear-wave splitting measurements were performed at nearby stations [14,15]. There are also independent indicators of the orientation of the maximum horizontal principal stress in southern California. The Central Valley and Oklahoma stations were chosen because they had relatively little sedimentary cover, and there are independent stress indicators available [14,16]. The Permian Basin Stations were chosen due to the economic significance of the region. In addition, the number of stations that we analyzed in each region was dictated largely by the computational costs of the inversions. These computational costs are related to the high quantity (>10,000) of receiver functions that were calculated with the total variation deconvolution (Appendix A). The disparity in the number of stations from the California regions to the stable intraplate regions is because we anticipated fewer tectonic impacts on the anisotropy in the latter regions, and thus tested more thoroughly.

The method that we used to assess the anisotropy of the crust using receiver functions is a nonlinear inversion to find a best-fitting model corresponding to the recorded receiver function data. This method has three distinct advantages over shear-wave splitting of Moho converted Ps phases, which is also found in the literature [10,17,18]. First, there is no need to subjectively isolate Ps phases, which despite quality controls, may introduce error into the shear-wave splitting method. Second, events without clear Ps phases can be used, which may contain valuable information about other time periods of the receiver function but may also complicate the stacking. The third advantage of this method is that it could in principle isolate the anisotropy to just the brittle, upper crust where the stress can be controlling the anisotropy.

The method we employed was developed originally by [19] and furthered by [20,21]. The first step in the methodology was to deconvolve the seismograms using the total variation method outlined in Appendix A. The azimuthal direction from which an earthquake arrives at the station (back azimuth) is important in the treatment of anisotropy because patterns in the amplitudes and arrival times of waveforms with different back azimuths can be used to calculate the velocity anisotropy beneath the station. Thus, we binned the receiver functions into 20-degree back azimuth bins and stacked the receiver functions within each bin after first correcting for the angle of incidence. The 20-degree bin is subjective and dependent upon the number of events per station. If the bin is too small, there are not enough events to stack. If it is too large, the sensitivity to back azimuth and thus anisotropy is lessened.

Once the raw data were processed for each station, the Ps phase was isolated, and all of the analyses presented below corresponded to inverting for just a narrow window of 2 s before and 2 s after the theoretical Ps arrival, which was acquired from the Earthscope Automated Receiver Survey for each location. The 2 s window was intended to capture the full Ps phase, which was anticipated to be most sensitive to crustal anisotropy. Next, the Neighborhood Algorithm was performed in order to find a best-fitting model for this nonlinear inversion [12,13]. The idea behind the Neighborhood Algorithm is that a set of forward models are calculated from an initial batch of models with randomly selected parameters. Next, a misfit calculation was performed for those initial forward models, and a Voronoi cell was placed around each model in parameter space such that when the next step was applied—a random walk from each of the initial models—the new models fell within the Voronoi cell. The misfit calculation was performed on all the new models. Then, the x number of models with the best misfits performed a random walk within their new Voronoi cells to create the next iteration of models. The parameter x was chosen by the user and was a function of both the desired completeness of search of the parameter space and the number of free parameters. The process iterated until convergence, which given the complexity of our model took about 40,000 models and 500 iterations for each station.

The forward modeling was performed using the software RaySum [22], which uses a ray theoretical algorithm to calculate synthetic seismograms for various back azimuths, dipping layers, and anisotropy. The most computationally expensive part of the inversion was performing the total variation deconvolution on the many thousands of seismograms generated by RaySum (Appendix A). The misfit calculation was performed by calculating a zero-lag cross-correlation coefficient between each modeled receiver function and the data with the corresponding back azimuth range within the four second Ps window. We used this coefficient averaged over every back azimuth bin to assess the goodness of fit, but because the Neighborhood Algorithm minimizes misfits, we had to subtract that value from one in the inversion.

The model we used for the forward calculations has three layers: the mantle, the lower crust, and the upper crust (Figure 2). The motivation for separating the crust into the lower and the upper is that the stress-controlled anisotropy is the result of the closure and dilation of fractures, which would only occur in the upper, brittle portion of the crust. For our objectives it was necessary to split the crust in this fashion. We inverted for a model with only seven parameters (

Table 1. Ranges of values searched for each free parameter.

Parameter	Upper Crust		Lower Crust		Mantle
Anisotropy Trend	0–180	Degree	-		-
Anisotropy (S)	0–10	%	0–20	%
Shear velocity Vs	3000–4000	m/s	3000–4000	m/s	4200–5200	m/s
Lower Crustal Anisotropy Trend	-		0–360	degree	-
Lower Crustal Anisotropy Plunge	-		0–90	degree	-

) by prescribing values for the thicknesses of the upper and lower crust obtained from the Earthscope Automated Receiver Survey H-k stacking [23]. We used a fixed relation between the P- and S-wave anisotropy ($1.64 An{i}_P=An{i}_S)$ and an empirical formula for the density of each of the free layers ($\rho =1.6612V_P-0.4721V_P^2+0.0671V_P^3-0.0043V_P^4+0.000106V_P^5)$ [24], and by assuming the $V_P=\sqrt{3}{V}_S\approx 1.73V_S$, which is the case for a Poisson solid. For a detailed argument on why fixing the P-wave anisotropy to the S-wave anisotropy is valid, see Appendix B. In this design, we assumed averaged mechanical properties for the entire crust, and the brittle/ductile transition mattered only for the anisotropy. Although eighteen parameters were used in the forward modeling, only the fast polarization direction of the upper crust (ϕ), the anisotropy percentage of the S-waves of each of the lower and upper crust (α_s), the trend and plunge of the anisotropy of the lower crust, and the shear velocity of the averaged crust and uppermost mantle (V_s) were inverted for, because the other parameters were either calculated from those seven or obtained from a reference. The Neighborhood Algorithm only searches over user-prescribed ranges of values for each free parameter. Those ranges are given in Table 1. The output of each inversion is tens of thousands of models with corresponding misfits. These were ordered by ascending misfit and the lowest-misfit model was assessed for the brittle fast polarization direction.

Testing the Inverse Algorithm

We tested the inverse algorithm using the deconvolution scheme for the receiver functions described in Appendix A against data from a synthetic model. For the purposes of the synthetic test, we used a model with anisotropy in both the upper and the lower crust, such that the successful results lend credence to the ability to differentiate between anisotropy in the two layers. The inputs and inverted parameters are shown in

Table 2. Results from the synthetic test.

	Average Cross-Correlation Coefficient [Unitless]	Upper Crustal Fast Polarization Direction [Degrees]	%S Anisotropy Upper Crustal [Percentage]	Vs Crustal [m/s]	Vs Mantle [m/s]	%S Anisotropy Lower Crust [Percentage]	Anisotropy Trend Lower Crust [Degree]	Anisotropy Plunge Lower Crust [Degree]
Synthetic Model	n/a	65	4.0	3222	4915	16.5	117	76
Inverted Model	0.97	64	4.0	3205	5192	17.2	116	76

. We chose to plot the model parameters in pairs of similar physical significance (e.g., V_s crustal and V_s mantle), color-coded by the lowest-misfit model in each region of parameter space (Figure 3). These plots allow for the quick assessment of where in the user-prescribed parameter space the best models reside. That is, in the dark squares in Figure 3, the misfit is low and there are good models. In the bright green squares, the misfit is high, and the Neighborhood Algorithm will avoid these regions. The four plots corresponding to the seven parameters can be seen in Figure 3 (the % S-wave anisotropy is plotted twice for visualization purposes—there are an odd number of parameters). The inversion captures all the parameters, and the waveforms were recovered effectively (Figure 4). In nearly every parameter, an almost perfect fit was observed between the best-fitting model and the known input. Thus, for this test, the inversions performed effectively.

3. Results

It is instructive to assess the signal quality of the processed data. There is often poor signal-to-noise ratio in the data, so the results for most of these stations could not be used (

Table 3. The qualitative signal-to-noise ratio and zero-lag cross-correlation coefficient for each station used in this study. The × marks indicate stations with insufficient signal quality.

Region	Station	Average Zero-Lag Cross-Correlation Coefficient	Assessment of Signal Quality
Central Valley, CA	S05C	0.51	√
	T05C	-	×
	T06C	-	×
	U05C	-	×
	PKD	0.74	√
Oklahoma	T33A	0.68	√
	T34A	0.72	√
	U33A	-	×
	U34A U35A V32A V33A V34A V35A W32A W33A W34A X32A X33A X35A	0.71 0.71 - - - 0.66 0.81 - - - - 0.72	√ √ × × × √ √ × × × × √
Southern California	DEC	-	×
	BFS	0.63	√
	BBR	0.74	√
	PLM	0.81	√
Permian Basin, TX	128A	-	×
	129A 130A 131A 227A 228A	- 0.72 - 0.49 0.48	× √ × √ √
	229A	0.65	√
	230A 231A 326A 327A 328A 329A 330A 331A 428A 429A	- - - 0.52 - 0.48 0.60 - - 0.55	× × × √ × √ √ × × √

). Therefore, we both qualitatively assessed the signal-to-noise ratio of the stacked data (Figure 5), differentiated by back-azimuth, for each station, and measured the average zero-lag cross-correlation coefficient for each of the best-fitting models relative to the data. The qualitative assessment of signal-to-noise ratio was quality control, in which we checked for the presence of distinct P and Ps phases in the receiver functions. We discarded any results where there were not distinct P and Ps phases or there was an average cross-correlation coefficient value below 0.30.

The average zero-lag cross-correlation coefficients of Table 3 are an indication of the quality of the inversion, yet it is helpful as well to visualize the results. Figure 6 shows the raw data as well as the waveforms from the best-fitting model for the PLM station in southern California. According to the cross-correlation coefficient, this is the highest quality of all of the twenty inversions, although station W32A in Oklahoma has nearly the same cross-correlation coefficient. Visual inspection also corroborates that the inversion was successful. While the waveforms from the data are not captured perfectly, the match is sufficient to validate the efficacy of our methodology.

Table 4. Full results from the best-fitting model for each of the twenty stations.

	Station	Fast Pol Dir	% Upper	Vs Crustal	Vs Mantle	% Lower	Trend	Plunge	Ave-Cross-Corr
Central Valley	PKD	159.54	2.18	3040.97	4472.89	9.61	347.20	83.28	0.74
	S05C	107.27	9.94	3000.08	4371.23	7.43	200.61	22.52	0.52
Oklahoma	T33A	118.50	9.98	3001.45	4828.87	9.11	15.78	65.60	0.68
	T34A	104.07	8.26	3412.47	4842.84	16.75	238.28	49.24	0.72
	U34A	176.66	0.05	3186.66	4629.48	14.04	97.76	88.50	0.71
	U35A	111.94	9.96	3384.27	4278.47	8.08	36.73	14.15	0.71
	V35A	154.90	7.46	3877.31	4718.13	19.43	87.80	2.24	0.66
	W32A	41.57	9.99	3251.72	4301.47	7.77	37.33	11.62	0.81
	X35A	141.06	9.96	3109.19	4914.64	10.74	141.25	78.14	0.72
Permian Basin	130A	104.56	9.99	3197.42	4855.14	10.44	43.87	61.75	0.72
	227A	73.67	6.32	3000.45	5134.44	7.95	137.68	30.88	0.49
	228A	99.32	2.85	3000.39	5117.56	5.34	277.19	61.70	0.48
	229A	79.15	0.34	3711.20	4306.89	13.37	60.85	19.23	0.65
	327A	166.59	0.13	3256.53	4217.36	19.76	113.42	48.72	0.52
	330A	170.92	4.97	3573.39	4794.87	19.99	271.13	21.01	0.60
	429A	12.22	3.46	3368.83	5141.99	18.75	219.87	31.95	0.55
Southern California	BBR	173.63	0.58	3967.05	4646.42	11.49	66.11	86.46	0.74
	BFS	125.60	10.00	3990.80	4896.44	17.77	79.87	61.96	0.63
	PLM	130.38	9.01	3983.57	4434.02	8.91	40.00	55.89	0.81

shows the parameters of the best-fitting model for all twenty stations tested.

The results detail a comparison of the direction of the maximum horizontal principal stress, as understood from previous studies, to the fast polarization directions obtained from our inversions. The previous studies on maximum horizontal principal stress directions rely on several methods. First, there are earthquake focal plane mechanism inversions, which take regional data on the radiation patterns from earthquakes to invert for the state of stress that best describes the data [3]. Second, wellbore data, whether drilling-induced tensile fractures or wellbore breakouts, are used to provide this measurement [3]. Third, the alignment of hydraulic fractures, typically determined from microseismicity, is used [3]. Finally, shear-wave splitting measurements from local earthquakes constrain the direction of crustal anisotropy [14,15]. These measurements from previous studies are referred to here as independent measurements of S_Hmax.

We hypothesize that there are four separate phenomena represented by the fast polarization directions of our inversion. First, the stresses are acting on quasi-randomly oriented fractures and faults, inducing a stress-controlled shear-wave velocity anisotropy. Second, there are some stations that lie close to large, crustal-scale discontinuities such as the San Andreas or Elsinore faults, and as was shown in [9], these features tend to align the shear-wave velocity anisotropy to the strike of the discontinuity. Third, the dip of the Moho can often cause effects on receiver functions that are nearly indistinguishable from anisotropy [25]. Finally, the inversion has been observed to be less accurate at high values of lower crustal anisotropy (Appendix C). Therefore, we see the combination of the first three factors in each of the fast polarization directions, and we see the fourth factor in the subset of the polarization directions that have high lower crustal anisotropy.

In Figure 7, we see the effect of the data quality wherein only two out of the five stations were able to be analyzed. Indeed, the S05C fast polarization direction aligns poorly with the independent measurement of S_Hmas, as noted by [3]. In contrast, the data from the PKD station along the San Andreas fault do not align to the known independent stress indicators, which means that that anisotropy is not stress controlled. The PKD inversion result does align with the strike of the San Andreas Fault, however, which may mean that the inversion has identified structural anisotropy. The blue arrow indicates an interpretation of the dip of the Moho [26]. The reason for the small number of stations tested (five) is due to the computational cost of these inversions, and we decided to use the majority of our resources in Texas and Oklahoma where there is not expected to be strong structural anisotropy.

In Figure 8, there are three fast polarization directions that disagree not only with the regional independent measurements of S_Hmax but also indicators of the regional upper crustal anisotropy [14,15]. This means that the inversion is not obtaining the upper crustal anisotropy for this region, which could be the result of the other three factors influencing the inversion: structural anisotropy, Moho dip, and the limitations of the inversion at high degrees of lower crustal anisotropy. The blue arrow indicates an interpretation of the Moho dip from [27,28], and there is a correlation between that direction of dip and the anisotropy indicators of stations BFS and BBR. The PLM station overlies the Elsinore Fault and aligns well to the strike, which could be an indication of structural anisotropy. Once again, computational cost was the primary impediment to analyzing more stations.

In Figure 9, there are eight out of seventeen stations that had both a qualitatively significant signal-to-noise ratio, and an average zero-lag cross-correlation coefficient greater than 0.30. The average of those coefficients is 0.56. Five of the seven fast polarization directions agree well with the known stress field. The other two (330A and 429A) seem to align at least partially with the Moho dip, and furthermore have lower crustal anisotropy values greater than 18% (see Table 4 for the full best-fitting model results). This mixture of Moho dip, high lower crustal anisotropy, and stress alignment may indicate that in certain subregions the effects from the upper crustal anisotropy might be more pronounced and lead to the good agreement with the stress field, while in other subregions the other three factors are dominant. The Moho dip is interpreted for the Permian basin from Shen and Ritzwoller (2016) [29].

Figure 10 shows some of the complexities of using teleseismic waves to determine the upper crustal stress orientation. Of the fifteen stations analyzed here, seven have sufficient signal. On average, the average zero-lag cross-correlation coefficient is 0.72, which is considerably higher than that of the Permian basin, most likely due to reverberations in the Permian basin sedimentary layers decreasing the signal-to-noise ratio of the receiver functions. Three of the stations in the northern portion of the study region (T33A, T34A, and U35A) exhibit a fast polarization direction that is relatively close to the independent measurements of S_Hmax. The degree of lower crustal anisotropy according to the best-fitting model for the V35A station is significant (>19%), but the rest of the stations do not have extreme values for the lower crustal anisotropy. We posit that either complicated structural anisotropy due to faulting in the region or fabrics in the lower crust that have been uplifted into the upper crust are the most likely explanations for the moderate performance of the inversion. The Moho dip is not strong or systematic in this region and has thus been omitted from the map [29].

4. Discussion

Refs. [19,21] successfully apply techniques similar to the nonlinear inversion method of this paper to regions in the Tibetan Plateau and Parkfield, California, respectively. In both studies, the authors were interested in anisotropy of the lower crust. Indeed, the potential degree of anisotropy due to mineral alignment may be as high as 20% [11], which is significantly greater than the 2–4% reported for stress-controlled anisotropy in the upper crust [8,9]. Therefore, the upper crustal signature that we want to discern is smaller than the lower crustal and mantle anisotropy wherein there is a greater potential for mineral alignment due to ductile flow. Even though the nonlinear inversion allows us to search for just the brittle crustal anisotropy, the significant contribution from the lower crust still had the potential to obfuscate the results (Appendix C). This is evidenced in our results by the fact that our measurement at the BFS station (Figure 8) did not align to the anisotropy of [14,15]. From this result, it seems clear that the inversion only sometimes successfully differentiated anisotropy between the upper and the lower crust. The other factor that appears to play a role in the results from the nonlinear inversion is the direction and magnitude of the Moho dip. In Southern California (Figure 8), the BFS and BBR stations both seem to align to a steeply dipping Moho, while in the Central Valley (Figure 7) there is no alignment to the Moho dip at roughly 3°. In the Permian Basin, where the dip magnitude is only roughly 1.5°, there are only potentially two stations influenced by the Moho dip (330A and 429A). We must emphasize, however, that the Moho dip depicted in Figure 7, Figure 8 and Figure 9 is only an interpretation of published data, and the true nature of the dip is likely more complicated.

Even in regions where the interplay of lower crustal anisotropy may be too significant to deduce a reliable upper crustal signal, as is seen in many of the measurements of this study, there is still information about the crust contained in these inversions. The measurements that do not align to the direction of the maximum horizontal principal stress may indicate that the structure of the crust is complicated and may be modeled with more sophisticated methods that rely on plunging-axis anisotropy [30]. The results that differ may also represent locations where fabric from the lower crust has been uplifted into the upper crust. The validity of the results that align to the direction of the maximum horizontal principal stress is enhanced by those results that indicate a more complicated structure, in that both reveal different, yet interrelated, properties of the crust. A challenge of this technique is that it may not be possible to discern between complicated structure and stress without prior knowledge of the stress state or crustal structure.

One of the potential problems that our algorithm faces is that the anisotropy is only allowed to have a tilted axis in the lower crust and not the upper crust. The reason for this is that we are interested in the component of anisotropy that lies in the horizontal plane and may be the result of tectonic stresses. Furthermore, it has been shown that often the strongest signal in a receiver function is the two-lobed pattern with back azimuth as a result of tilted-axis anisotropy [31,32,33], yet it is unclear whether our synthetic receiver functions exhibit this behavior. In addition, our model has significant constraints placed on it, such as the two-layered system in the crust, which cannot incorporate multiple, thin anisotropic layers [32]. These constraints also limit the ability to analyze fabrics developed in the lower crust that are later uplifted into the upper crust. The result of the incomplete coverage of tilted-axis anisotropy as well as the constraints on the model means that our algorithm, while it likely resolves a portion of stress-induced anisotropy, does not address the full range of realistic crustal anisotropy.

There are several high values of crustal shear-wave velocity in southern California from our inversion results (Table 4). These may be the result of the restrictions on the model. It is also possible that the difference in the values correspond to different depths sampled by the modeled Ps phase. In a previous study in southern California, with a depth of investigation of 10 km, similar values to our inversion were seen [34]. Furthermore, high S-wave velocities are seen in locations throughout the North American crust [35].

There is a discrepancy in the average zero-lag cross-correlation coefficient values between the Permian Basin and Oklahoma (Table 3). That is, the Permian Basin values are on average smaller than those of Oklahoma. This is likely evidence for basin multiples interfering with the receiver functions in the Permian Basin, but the relatively less sedimentary cover in Oklahoma is allowing for a cleaner, lower misfit inversion. Therefore, when using the nonlinear inversion method, basin multiples could play an important role in the quality of the results.

5. Conclusions

The interplay between the upper and lower crustal anisotropy proved to cause a clean, solely upper crustal signal to be difficult to discern. The upper crustal signal in the nonlinear inversion method was obfuscated by lower crustal anisotropy, as is evidenced by the misalignment between the inversion results and the previously published upper crustal anisotropy from shear-wave splitting measurements on local earthquakes in southern California. In the Permian Basin, five out of the seven stations that could be tested aligned to the stress field, which may indicate that the structural complexity and the interference from the lower crustal anisotropy is less than, for instance, in Oklahoma, where the majority of the stations show upper crustal anisotropy that is misaligned to the stress field. The data quality was an issue in this study, where only twenty of the forty-one stations analyzed had both a good signal quality and a high average zero-lag cross-correlation coefficient. Those stations that did indeed align to the independent measurements of S_Hmax lend credence that perhaps if the lower crustal anisotropy is weak, and the upper crustal anisotropy is both appreciable and the result of stress, the nonlinear inversion method could be used to determine the direction of S_Hmax. Finally, a new method of receiver function deconvolution was developed, which creates more impulsive receiver functions, but whose impact on this study was difficult to discern.