A 3D Full Stress Tensor Model for Oklahoma

The stress tensor is an important property for upper crustal studies such as those that involve pore fluids and earthquake hazards. At tectonic plate scale, plate boundary forces and mantle convection are the primary drivers of the stress field. In many local settings (10–100 s of km and <10 km depth) in tectonic plate interiors, we can simplify by assuming a constant background stress field that is perturbed by local heterogeneity in density and elasticity. Local stress orientation and sometimes magnitude can be estimated from earthquake and borehole‐based observations when available. Modeling of the local stress field often involves interpolating sparse observations. We present a new method to estimate the 3D stress field in the upper crust and demonstrate it for Oklahoma. We created a 3D material model by inverting multiple types of geophysical observations simultaneously. Integrating surface‐wave dispersion, local travel times and gravity observations produces a model of P‐wave velocity, S‐wave velocity, and density. The stress field can then be modeled using finite element simulations. The simulations are performed using our simplified view of the local stress field as the sum of a constant background stress field that is perturbed by local density and elasticity heterogeneity and gravitational body forces. An orientation of N82°E, for the maximum compressive tectonic force, best agrees with previously observed stress orientations and faulting types in Oklahoma. The gravitational contribution of the horizontal stress field has a magnitude comparable to the tectonic contribution for the upper 5 km of the subsurface.

• We compute the stress field using finite element modeling by considering contributions from both gravitational and far-field tectonic forces • We obtain the orientation and magnitude of the tectonic force by finding the stress field that best-fits the observed stress observations • Gravitational contribution to the horizontal stress field has a comparable magnitude regarding tectonic contribution for the upper 5 km Supporting Information: • Supporting Information S1 • Data Set S1 • Data Set S2 • Movie S1 • Movie S2 • Movie S3 Correspondence to: C. Chai and A. A. Delorey, chaic@ornl.gov; andrew.delorey@lanl.gov Heidbach, 2014) have been used to explain large-scale deformation (e.g., mountain building) and seismicity (both interplate and intraplate). The importance of the gravitational contribution in stress modeling has been emphasized by Humphreys and Coblentz (2007) and Flesch et al. (2000Flesch et al. ( , 2007. Forte et al. (2007), Levandowski et al. (2016) and Lund Snee & Zoback (2020) suggest both tectonic forces and subsurface elastic property heterogeneity may have contributed to the seismicity of the New Madrid seismic zone. Flesch et al. (2007) found that the gravitational potential energy and tectonic forces contribute to the deviatoric stress magnitude almost equally in western North America. However, it is not clear whether the gravitational and tectonic contributions are similar to the deviatoric stress magnitude in other regions. We chose the Oklahoma area to study the 3D stress field variations since previously published stress observations in the area can validate our stress modeling method. We hypothesize that most of the observed deviations from the tectonic stress field in Oklahoma are due to gravitational effects and heterogeneity in material properties including density and elastic modulus. By modeling both the gravitational and tectonic components of the stress field, we can both test this hypothesis and constrain principal component magnitudes. For the purposes of this paper, "tectonic" stress refers to stresses that are imparted on our study area through boundary conditions and includes such things as ridge push, slab pull, and basal traction. More detailed material models have been suggested as one approach to improve stress models by several researchers (e.g., Ghosh et al., 2019;Reiter & Heidbach, 2014). In order to model the gravitational and tectonic components of the stress field in Oklahoma accurately, we require a detailed 3D material model and information (direction and magnitude) about tectonic forces. Since directions and magnitude of tectonic forces are not known for most regions, we estimate them with stress observations via multiple stress modeling tests. The material model can be obtained from subsurface seismic structure investigations.
A few detailed investigations have focused on the 3D crustal structure in the Oklahoma region. Evanzia et al. (2014) studied the upper mantle structure using teleseismic body wave tomography in the Texas and Oklahoma area. A few regional and continental scale studies mainly based on Earthscope Transportable Array data with nominal station spacing of ∼70 km (e.g., Chai et al., 2015;Porter et al., 2016;Schmandt et al., 2015;Shen & Ritzwoller, 2016) imaged the broad-scale heterogeneity of the Oklahoma region. Limited by different scope and focus of these studies, fine scale subsurface structure of both P-and S-wave velocities in the study area is not well resolved. Seismic velocity models produced by these broader scale studies can be used as a starting model for more detailed investigations. The compiled global model Crust 1.0 (Laske et al., 2013) was produced by combining Earth models from active source surveys. However, the coverage and resolution of these earlier studies are not sufficient for our stress modeling. We constructed a material model to extend these earlier efforts.
In this study, we first present data used in our technique that includes observations for material model inversion and for stress modeling. Details on the inversion for the material model and stress modeling are documented in the following section. Then, we present results for the material model and the stress model of the Oklahoma region. The major findings, assumptions, and limitations are discussed in the last section.

Data
Our analysis consists of constraining a material model (P-and S-wave velocity and density) and stress modeling, which requires both geophysical observations and stress measurements.

Observations for Material Model Inversion
We simultaneously and jointly invert local P-and S-wave first arrival times, Rayleigh-and Love-wave dispersion maps, and satellite-derived Bouguer gravity observations to model the 3D subsurface P-and S-wave velocity and density structure. Each of these data sets went through various quality control procedures.
Interactive visualization tools (Chai et al., 2018) were used to visually examine the datasets as well as the resulting material models. Details of the data selection and preprocessing are documented in the following paragraphs.
The local body waves arrival time data set (catalog locations and P-and S-arrivals) was obtained from United States Geological Survey (McNamara et al., 2015). Catalog arrival times that we suspected to be erroneous, such as those leading to more than 5 s misfits (unrealistically large), were excluded from the inversion for material properties. Only events with six or more P-wave arrival times are included in the inversion so the location of seismic events is well constrained. The final arrival times data set used in the inversion comprises 58,896 P-wave arrivals and 20,155 S-wave arrivals from 3,740 earthquakes that were recorded at 157 seismic stations between November 2011 and July 2015 ( Figure 1a). Besides absolute catalog body wave arrival times, differential catalog times were computed (Waldhauser, 2000) and included in the inversion (Zhang & Thurber, 2003) because of their ability to further improve the relative locations of neighboring earthquakes and their surrounding velocities. As shown in Figure 1, most of the earthquakes are located in north central Oklahoma, and the station distribution is reasonable with a good azimuthal coverage.
The surface wave dispersion estimates (used as observations in this study) were obtained from an update of Herrmann et al. (2016). The dispersion model was constrained with local earthquake and ambient noise measurements. Group and phase velocities of both Rayleigh and Love waves are available online (Herrmann et al., 2016). We use periods ranging from 5.0 to 51.5 s with 24 data points per dispersion curve when sufficient data are available. The surface wave data were resampled from the dispersion model onto a 0.5° by 0.5° horizontal grid that covers the study area. Four dispersion curves at each of these horizontal grid nodes are included in the inversion when available. The Bouguer gravity data (observations for this study) was extracted from the global model WGM2012 (Balmino et al., 2012) and sampled to the same horizontal grid as for the dispersion observations. The gravity data was filtered in the wavelength domain to avoid overfitting and to remove long-wavelength features mainly caused by deep density variations. Very short wavelength (less than 50 km) and long wavelength (longer than 200 km) anomalies were excluded since they are beyond our spatial resolution (50 km laterally). Gravity observations at the four edges (around 200 km wide) were not used in the inversion to avoid edge effects (no gravity contributions outside of the study area). We used a larger area (the entire region as shown in Figure 1a) for the velocity model inversion so that the edge effects do not influence our results.

Observations for Stress Modeling
During the stress modeling, we used S Hmax orientations and A ϕ (Simpson, 1997) as data constraints. A total of 94 S Hmax orientations with uncertainties of 15°, 20° or 25° referred respectively as qualities A, B, and C (Heidbach et al., 2018) were accepted as stress observations (see Figure S1). A ϕ (Simpson, 1997) is determined from the magnitudes of the three principal stresses and ranges from 0 (radial extension) to 3 CHAI ET AL.

10.1029/2020JB021113
3 of 14 (radial constriction), where the values for pure normal, pure strike-slip, and pure reverse are 0.5, 1.5, and 2.5, respectively. We used A ϕ values from Levandowski et al. (2018) in our stress inversion because they cover our study area and have published uncertainties, though other interpretations exist (Lung Snee & Zoback, 2016).

Method
We first describe the approach that we used to obtain the material model. Then, we document details on the stress modeling.

Inversion for the Material Model
The joint inversion technique we used is the same as that of Syracuse et al. (2016Syracuse et al. ( , 2017 which used multiple geophysical datasets to directly constrain P-and S-wave velocities, and mass density. The inversion algorithm was developed based on tomoDD (Zhang & Thurber, 2003, 2006 and an inversion program that simultaneously inverts surface-wave and gravity data (Maceira & Ammon, 2009). The shear velocity-density relationship from Maceira and Ammon (2009) was used to relate gravity observations to seismic speeds. We started the inversion using a 3D initial model from Chai (2017), which was constrained with spatially smoothed P-wave receiver functions (Chai et al., 2015), surface-wave dispersion (Herrmann et al., 2016) and gravity observations (Balmino et al., 2012). The inversion iteratively updated the P-and S-wave velocity models by minimizing the data fits and taking into consideration other regularizations constraints such as smoothing and damping. Different weights were assigned to each type of data constraints. We performed two suites of inversions to find the optimal combination of weights using L-curve analysis as suggested by Syracuse et al. (2015). The first set of inversions searched for the appropriate damping and smoothing values (see Figure S2 and Visualization S1 in the electronic supplements). The second set explored the weights for gravity and surface-wave data (see Figure S3 and Visualization S2 in the electronic supplements). The optimal material model (P-and S-wave velocities and density) was selected by minimizing an objective function of data fits and regularization terms similar to Syracuse et al. (2016Syracuse et al. ( , 2017. The optimal material model fitted the travel times (Figure 2), gravity observations (Figure 3), and surface-wave dispersions (Figures 4 and 5 and S4) reasonably well. Earthquake locations were updated during the inversion simultaneously as well.

Stress Modeling
We used a commercial finite element code (ABAQUS standard/implicit, version 2018) to model both the tectonic and gravitational body force components of the stress field. The full finite element model is 2,000 km along each horizontal dimension and 200 km in depth (Figure 1b). The model was designed to accurately represent stress in the uppermost crust (<10 km), and we do not interpret any results for the deeper crust or upper mantle. Our study area, which falls geographically within 33.5-37.0 N latitude and 100 to 94.5 W longitude and is approximately 390 km south to north and 500 km west to east, is centrally embedded within the full model. The axes of this central part of the model were oriented west-east and south-north while the full model was rotated by 0°, 8°, 12°, 16°, 20°, or 24° counterclockwise to correspond to the principal orientations of the modeled tectonic stress field. These directions correspond to azimuths of N90°E, N82°E, N78°E, N74°E, N70°E, and N66°E, respectively, mimicking the general ENE-WSW S Hmax orientations in central and eastern North America (Alt & Zoback, 2017;Heidbach et al., 2018;Levandowski et al., 2018;Walsh & Zoback, 2016). The source of this regional stress field is not the focus of this study, but such directions closely accord with either Mid-Atlantic Ridge push or basal drag of North American absolute motion (M. L. Zoback & Zoback, 1980).
The finite elements are 8 node hex, elastic elements with horizontal dimensions of 12 km by 12 km and vertical dimensions of 2 km near the surface, increasing to 20 km vertical spacing near the base ( Figure S5). The material properties (P-and S-wave velocities, and density) from the joint inversion were mapped to the geographically appropriate elements. Elements that fall outside the boundaries of our material model were assigned the median value for their depth. Assuming the elastic model is CHAI ET AL.
10.1029/2020JB021113 5 of 14  isotropic, we converted the P-and S-wave velocities, and density into Young's modulus and Poisson's ratio (input for ABAQUS) using formulas from Shearer (2009). The formulas are also included in supporting information S1.
We modeled the effect of body forces (gravity) by holding constant normal displacement boundary conditions on the four sides and bottom and then apply the gravitational field. This process was performed using a "geostatic" step in ABAQUS. Normally, a volume with body forces applied and zero normal displacement boundary conditions on the bottom and all horizontal boundaries would compress under its own weight.
In this geostatic procedure, we determined equilibrium stress conditions for the model due to the applied body forces and boundary conditions, while maintaining the original dimensions. With this procedure, we began with a model that is pre-stressed with overburden before we apply tectonic stress. The overburden stress is approximately the product of density and gravitational acceleration integrated over the depth. Due to 3D heterogeneity in density and elasticity, the stress field caused by overburden is also heterogeneous.
CHAI ET AL.
10.1029/2020JB021113 6 of 14 We modeled tectonic stress by imposing a small, uniform (inward) displacement of 61 m along the southwest facing boundary of the full model and zero displacement boundary conditions on the northwest, southeast, and northeast sides (Figure 1b). The displacement applied is approximately equivalent to a tangential compressional stress of 1.9 MPa in the direction of the displacement in the upper 5 km. Since the elastic model is linear, we increment the resulting stress field to model different tectonic stress magnitudes. We intend the applied boundary conditions to approximate any regional contributions to the stress field including plate boundary forces and uniform tractions at the bottom of the lithosphere.
We determined the magnitude of the tectonic force for the best overall fit to the available S Hmax orientation observations and reported A ϕ values. The depth of measured S Hmax orientations varies, with those derived from borehole measurements at 1 km or more and earthquake moment tensor solutions at 4-5 km. To determine fit to S Hmax orientations, we compared the modeled S Hmax orientations with observed S Hmax orientations ( Figure 9). We also used faulting type as a constraint because faulting type is determined by the relative magnitudes of the principal stresses (e.g., Simpson, 1997). In our study area, the dominant style of faulting is strike-slip (Dziewonski et al., 1981), but there are significant and spatially variable amounts of net horizontal extension accommodated by oblique-normal and normal faulting (Levandowski et al., 2018;Walsh & Zoback, 2016).
Using the orientation of 94 available stress indicators (S Hmax ) and the faulting types in Oklahoma, we search for the orientation and magnitude of the tectonic force by finding the model that best fits (with the minimum misfit) the stress observations. When both the magnitude and orientation of S Hmax are available, alternative approaches (Reiter & Heidbach, 2014) can also be used to determine the orientation and magnitude of the tectonic force. We calculated all plausible magnitudes for a set of orientations and identified the best magnitude and orientation of the tectonic force by minimizing the following function where G is the Jacobian matrix describing the relationship between the model parameters and the model predictions, m is model vector, C d are the reported data uncertainties, and d is the data vector. The model vector has two values, the tectonic force magnitude and orientation. The data vector is composed of observed stress orientations (S Hmax ), and A ϕ . The misfit of each prediction is scaled by the uncertainty of the underlying observation using C d . S Hmax is a measure of the relative magnitudes of the two horizontal principal stresses, while A ϕ is a measure of the relative magnitudes of all three principal stresses. Combining these two types of stress observations provide constraints on both the direction and magnitude of the tectonic force.

Results
Our results include a 3D material model and a 3D stress model.

Material Model
We use depth slices, cross-sections, and interactive visualization to present the 3D material model. Depth slices and two cross-sections of the S-wave velocity model are shown in Figures 6 and 7, respectively. Depth slices of the P-wave velocity and density models are shown in Figures S6 and S7 in the electronic supplements, respectively. An interactive tool (Chai et al., 2018) is provided in the electronic supplements (Figure S8 and Visualization S3) to view the seismic velocities at other locations in the model. Our material model confirms previous geophysical and geological surveys but with a more complete image. At 3 km depth (Figure 6), the slow anomaly in southwest Oklahoma corresponds to the Anadarko basin. The anomaly extends more than 10 km in depth as we can see in the cross-section A1-A2 (Figure 7), which agrees with Johnson's geological compilation (Johnson, 2008). Low seismic velocities (upper 10 km in Figure 7) are likely associated with the Anadarko basin. The relative location between the imaged low velocities and the majority of the earthquakes is consistent with Isken and Mooney (2017)'s results. In the mid and lower crust (∼15-40 km in depth), we image a high-velocity anomaly beneath the Anadarko basin. The high-velocity anomaly may be related to the past igneous activity associated with the Southern Oklahoma Aulacogen (Gilbert, 1983). We found a slower upper mantle beneath the southern Oklahoma that is consistent with Evanzia et al. (2014)'s model.

Stress Model
The orientation of the observed S Hmax for the majority of the stress indicators varies between East and N60°E. In an attempt to constraint the orientation of the tectonic forces applied to our model, the misfit between the observed and predicted stress orientation was evaluated for a range of boundary force orientations which represent the far-field tectonic forces ( Figure 8). For all the cases considered, the misfit is high for very low magnitudes of tectonic force (displacements < 3 km). We interpret this high misfit as a consequence of gravitational body forces dominating the calculated stresses in the model. Conversely, the high misfit for very high tectonic forces (displacements >7 km) is the consequence of tectonic stresses dominating the predicted stress orientations. For most of the cases considered, the misfit is minimized for tectonic stresses imparted by a range of ∼4 to ∼6 km displacement. For orientations in the range of E-N66°E, the two best fitting realizations are models with a displacement of 3,965 m (∼123 MPa) with an orientation of N82°E and a displacement of 5,551 m (∼173 MPa) with an orientation of N70°E. Though the misfit is nearly identical in both cases, the relative contributions of the S Hmax and A ϕ to the overall misfit are different (Figure 8). We examine the spatial distribution of data match for both models to choose a preferred model (Figure 9). The N82°E model fits the spatial variations of the observed A ϕ a little better than the observed S Hmax directions (Figure 9). Compared to the N82ºE model, the N70°E model fits the observed S Hmax directions equally well, but not as well for the observed A ϕ (Figure 9). For this reason, we choose the N82ºE model as our preferred model. The fit to a few S Hmax directions may be improved with a material models at higher resolution.
Higher displacements produce models more toward reverse faulting and a more spatially uniform orientation for predicted S Hmax , and lower values CHAI ET AL.

10.1029/2020JB021113
8 of 14   Vs (km/s) produce a model more toward normal faulting and more variable spatial orientations for predicted S Hmax . An arbitrarily high tectonic force is not supported by the data both in terms of the orientation of S Hmax and the faulting type for the best fitting stress orientations (Figure 8). The best fitting model, that does not consider gravitational body forces, is for a stress orientation of N66°E with an arbitrary displacement magnitude, and has a normalized misfit of 3.85 compared to 1.01 for our preferred model (Equation 1). This means that stress caused by gravitational body forces explains most of the residual misfit of a uniform model, supporting our hypothesis that most of the observed deviations from the tectonic stress filed in Oklahoma are due to gravitational effects and heterogeneity in material properties including density and elastic modulus.
According to previous studies, A ϕ (Simpson, 1997) ranges from ∼1.0 in northernmost Oklahoma to 1.5 in central and southern Oklahoma (Levandowski et al., 2018;Lund Snee & Zoback, 2016;Walsh & Zoback, 2016). A ϕ value of 1.5 corresponds to strike-slip faulting. As A ϕ values increase from 1.5, the faulting type evolves toward reverse faulting and as values decrease from 1.5 faulting type evolves toward normal faulting. Though poorly constrained by focal mechanism inversions, portions of southern Oklahoma may have A ϕ of up to 1.7 (Levandowski et al., 2018;Lund Snee & Zoback, 2016). A ϕ of the N82°E stress model falls within the range for strike slip faulting with some regions in the southwest having higher values. The models with a tectonic force orientation of E and N66˚E have considerably higher misfits to the weighted sum of S Hmax and A ϕ (Figure 8). The resulting stress field has an S Hmax orientation from southwest to northeast, consistent with average S Hmax measurements (Alt & Zoback, 2017;Heidbach et al., 2018) and focal mechanism inversions (Levandowski et al., 2018;Walsh & Zoback, 2016). The deviatoric shear stress associated with these stress magnitudes are in the range of 10-30 MPa on optimally oriented faults at earthquake CHAI ET AL.
10.1029/2020JB021113 9 of 14 depths (around 5 km). The stress orientations and deviatoric stress magnitudes predicted here are consistent with the intraplate stress field predicted in Oklahoma (Alt & Zoback, 2017) and mid-plate North America as described in a number of previous studies; including qualitative assessments of the states of stress (e.g., Zoback, 1992;Zoback & Zoback, 1980), numerical models of the intraplate stress field using whole-plate linear elastic shell models (e.g., Humphreys & Coblentz, 2007;Richardson & Reding, 1991), as well as more recent evaluations of the relative stress magnitudes and orientations (Heidbach et al., 2010(Heidbach et al., , 2018(Heidbach et al., , 2007Lund Snee & Zoback, 2020). While the orientation of the S Hmax for the intraplate stress field is fairly well-constrained and understood in the context of local heterogeneities (see discussion in Schoenball & Davatzes, 2017), the magnitude of the intraplate stress field remains more elusive. The stress magnitudes predicted here are consistent with the body of evidence that the tectonic stresses in the lithosphere are generally of the order of tens of MPa averaged over the thickness of the lithosphere.
We examine the stress field at the 5 km depth for the Oklahoma area since most large earthquakes in the region occurred around that depth.
Comparing the horizontal component of the gravitational contribution with the tectonic contributions ( Figure 10), we found the gravitational contributions to the horizontal components of the deviatoric stress field (subtracting out the hydrostatic stress) is on the same order of magnitude to the tectonic contribution. The vertical component of the stress field is dominated by the gravitational contribution. Due to spatial changes in density and elastic properties, both the gravitational contributions and the tectonic contributions are non-uniform for the horizontal components. The spatial distributions of gravitational and tectonic contributions show different patterns. The resulting modeled deviatoric stress field shows spatial variability across Oklahoma. As a result, we need to include both the gravitational and tectonic contributions in the stress field calculation, which confirms previous studies for other regions (e.g., Levandowski et al., 2016;Reiter & Heidbach, 2014). The tectonic contribution to stress on optimally oriented faults is between 10 and 30 MPa within the study area; the gravitational contribution is 0-20 MPa, mostly of opposite sign, but of similar magnitude.
We acknowledge that tectonic forces acting on the study area through the boundaries are almost certainly not uniform in the Earth. But a uniform force at the boundaries is sufficient for testing the hypothesis that variations from a uniform stress field in the upper crust are primarily due to heterogeneous elasticity and density in the crust at local and regional scales. As our model is fully elastic, we are also neglecting the effect of viscoelastic behavior in the lower crust and upper mantle, which would have the effect of relaxing stresses at those depths with some time dependence. Our intent in this study is to present a method to model stresses in the upper crust where earthquakes and other activities occur, not calculate a stress model throughout the crust and upper mantle. Earth curvature will need to be considered when applying this method to a larger spatial scale. shows A ϕ (Simpson, 1997) values at 5 km depth computed from the stress model. A ϕ value of 1.5 corresponds to strike-slip faulting. As A ϕ values increase from 1.5, the faulting type evolves toward reverse faulting and as values decrease from 1.5 faulting type evolves toward normal faulting. The background in (a) shows A ϕ values (computed from measurements with a depth range of 0-16 km and an average depth of 5 km) from Levandowski et al. (2018).

Summary
We jointly inverted surface-wave dispersion, gravity, and local travel time observations for a 3D elastic property model for the Oklahoma region. The material model can be further improved with deep learning and transfer learning derived travel time observations (Chai et al., 2020). Utilizing the 3D material model, a model of the 3D stress tensor field for Oklahoma was computed by considering both gravitational and tectonic contributions. A model that includes both gravitational and tectonic components of the stress field fits observed stress indicators better than the tectonic component alone, indicating that the gravitational component helps to explain small-scale variations in principal stress orientations. We used observed stress indicators and faulting types to constrain the tectonic force orientation and magnitude equivalent. Our preferred model has a tectonic force orientation of N82°E and explains well the stress observations and the faulting types. An equivalent stress magnitude near 123 MPa (shortening of 3,965 m) fits the stress observations and the faulting types better than other magnitudes. This corresponds to a deviatoric shear stress of 10-30 MPa on optimally oriented faults. The stress field in the upper 5 km due to gravitational body forces has a comparable magnitude to the tectonic-driven stress field and the modeled stress field in the Oklahoma region has significant spatial variations.
CHAI ET AL.
10.1029/2020JB021113 11 of 14 Figure 10. A comparison of magnitudes of deviatoric stress (left column) and its gravitational (middle column) and tectonic contributions (right column) at 5 km depth using the preferred stress model. The color scale changes for each image. The deviatoric stress tensor D is the sum of the gravitational contribution G and the tectonic contribution T. The first subscript is the orientation of the surface that the tensor component acts on. The second subscript of the tensor component is the direction of component. Subscript 1 represents north, 2 for east, and 3 for down. Our results demonstrate that a reliable 3D stress field for the upper crust can be computed using a 3D material model and stress observations. When previous localized stress orientation measurements and focal mechanisms are known, the orientation and magnitude of the regional tectonic forces can be constrained through stress modeling. 3D stress modeling has significant advantages over traditional methods (e.g., interpolation, extrapolation) including the ability to obtain stress magnitudes and the ability to use measurements away from boreholes and earthquake focal mechanisms. Careful validations of 3D stress models are needed to ensure that the stress field is correctly simulated. As small-scale stress variations may favor or oppose future deformation, our technique can also help with subsurface engineering problems (with a smaller grid) and natural hazards. Though our model closely recovers most available stress indicators, the stress tensor continues to be a difficult property to validate, which could be done with new observations.