Teleseismic waves reveal anisotropic poroelastic response of wastewater disposal reservoir

1.   Introduction

Faults in areas of injection-induced seismicity, like central Oklahoma, have elevated pore pressure and as a result can have enhanced sensitivity to deformation. For instance an enhanced sensitivity to dynamic triggering in Oklahoma (van der Elst et al., 2013) may arise because pore fluid pressures are higher than they would be in a natural state of sub-hydrostatic equilibrium (Figure 1). In this case, and assuming differential stress and friction remain constant, the associated decrease in effective mean stress implies that relatively small stress perturbations — tens to hundreds of kPa — may be enough to cause frictional failure (e.g., Townend and Zoback, 2000). Dynamic stresses from large teleseisms can be on the order of this threshold. Current data on dynamically triggered earthquakes in areas of induced seismicity either cannot detect instantaneous triggering or don't look closely enough to observe it, which makes it difficult to say if triggering is due to shear or volumetric strain from Love or Rayleigh waves, respectively.

Figure  1.  Measurements of fluid pressure in the Arbuckle wastewater disposal reservoir in Oklahoma. (a) Formation fluid pressures in Osage County over time, relative to normal hydrostatic pressure. These data are based on drill stem tests (DSTs) (Hussey and Halihan, 2018) and continuous measurements in the Arbuckle (Barbour et al., 2019). For comparison we also show measurements at KGS well 1−28 in Sumner County, Kansas (Nolte, 2017). The inset map shows the locations of DSTs and continuous Arbuckle measurements in Osage County, along with seismicity and mapped faults from the Oklahoma Geological Survey (see Data and Resources). (b) Pressure measurements from DSTs in depth. For comparison we also show measurements from the Mervine Oil Field in Kay County, Oklahoma, the neighboring county to the west of Osage County (Davis, 1985). Recent measurements are closer to typical hydrostatic pressures (less sub-hydrostatic). For reference, the formation pressure at KGS 1−28 on 2016/07/30 was ~14.636 MPa; measurements are made at ~1539 m, where hydrostatic pressures should be ~15.098 MPa.

DownLoad: Full-Size Img PowerPoint

An additional deformation sensitivity that arises at elevated pore pressure is elastic anisotropy. Porous, fractured, and dilatant Earth materials such as porous sedimentary rocks, fault zones, and granular soils have elastic properties that depend non-linearly on the stress state at reduced effective stresses. For example stress-induced crack closure results in increasing Young's modulus and elastic wave speed with increasing effective confining pressure (Birch, 1960, 1961; Walsh, 1965a, b). The variation of moduli with effective stress requires that most Earth materials will have anisotropy whenever the differential stress is non-zero, especially when some form of damage (e.g., cracks) is introduced. That is, they are more compliant in the direction of the minimum principal stress than they are in the direction of the maximum principal stress, and will show sensitivity to changes in differential as well as mean stress. The specific response depends on the material, the stress state, and the sense and magnitude of the imposed changes in stress, the latter being the aspect to be exploited in the present study. At saturated undrained conditions the sensitivity is manifest as anisotropic poroelasticity.

In Southern California, dynamic triggering of microseismicity is common in regions of both natural and induced seismicity, but factors like peak ground velocity and dynamic strain are weak predictors of when triggering occurs (Fan et al., 2021). However, dynamic strains are useful for determining poroelastic response; isolating the response to shear and volumetric surface waves can provide constraints on material properties, pore pressure, and stress state. Thus, a richer understanding of how surface waves affect fluid-pressures in basal disposal reservoirs like the Arbuckle Group in Oklahoma could help understand the conditions for dynamic triggering, and how those are modified by fluid injection. For the purposes of simulating the effects of wastewater disposal, the subsurface is often simplified as a series of horizontal layers, each with isotropic poroelastic properties. In Oklahoma, it is unclear how appropriate the isotropic assumption is for the principal wastewater disposal reservoir — the Arbuckle reservoir — and whether anisotropic variations in material properties have a pronounced effect on the response to injection, and the generation of Coulomb failure stresses at seismogenic depths. In this paper we investigate the Arbuckle pore pressure response to dynamic shear and volumetric strains with surface waves from teleseisms.

2.   Methodology and Data Sources

Since 2017 the United States Geological Survey (USGS) has been monitoring fluid-levels in a repurposed Arbuckle disposal well in Osage County, Oklahoma (see insets in Figure 1b), with frequency sampling up to 1/4 Hz. Based on the tidal signals seen in this well, Wang CY et al. (2018) showed that leakage out of the reservoir is significant; Barbour et al. (2019) came to a similar conclusion, but they also showed a clear response to dynamic strains from the passage of long period surface waves, with peak dynamic fluid-levels appearing to be proportional to the magnitude and distance from the source (see their Equation (3)). In other words, rather than being perfectly confined by the layers above (e.g., shales and sandstones) and below (the crystalline basement), the Arbuckle reservoir shows evidence of being semi-confined, with fluid flow occurring towards the surface and/or into the basement.

Fluid-saturated reservoir rock that behaves as a linear, isotropic, porous medium should experience only pore fluid pressure changes from dilatational strains found in Rayleigh waves; it should be insensitive to the transverse shear strains found in Love waves. As Skempton (1954) showed, the change in pore fluid pressure in an undrained poroelastic material depends on mean stress and differential stress with constants $B$ and $A$, respectively:

where $\Delta \sigma_{\rm m} = (\Delta \sigma_1 + \Delta \sigma_2 + \Delta \sigma_3)/3$, $\Delta \sigma_{\rm d} = \left(\Delta \sigma_1 - \Delta \sigma_3\right)/2$, and $\sigma_1$, $\sigma_2$ and $\sigma_3$ are the principal stresses. In a linear poroelastic medium, $A = 1/3$ and there is no response to shear deformation; however, since Wang CY et al. (2009), numerous studies have documented seismically induced fluid pressure oscillations with an apparent response to shear strain waves. Here, we reanalyze the seismic records from Barbour et al. (2019) and collect new records, finding that the Arbuckle reservoir also shows a shear response. One clear example of this is from the 2018 M_w7.9 Kodiak Alaska earthquake (event id us2000cmy3) because the Love and Rayleigh waves are well-separated: the hydroseismic record shows signals that correspond to the S-wave and Love wave, in addition to the Rayleigh wave (Figure 2).

Figure  2.  Shear strain response of the Arbuckle reservoir in Oklahoma to the 2018 M_w7.9 Kodiak, Alaska, earthquake ~ 40° away. Left: azimuthal equidistant map of the epicenter centered on the Arbuckle well with plate boundary faults (dark lines). Right: timeseries of the fluctuations in fluid pressure in the reservoir compared to areal and shear strains inferred from radial and transverse ground velocities from the colocated broadband seismometer at the surface (see Methods).

DownLoad: Full-Size Img PowerPoint

Interpreting hydroseismic measurements from wells depends on an understanding of the flow regime in the reservoir and in proximity to the well, along with knowledge of the dynamic strain driving flow. This can be challenging because the response can be affected by the hydraulic and elastic properties of the reservoir rock. Permeability is a principal control on the hydraulic diffusivity and fluid flow in the crust, but it varies over orders of magnitude (e.g., Manning and Ingebritsen, 1999) and can change in time (e.g., Elkhoury et al., 2006). Natural and faulting-related damage, microcracks, and fracture or karst networks affect not only permeability but also the elastic properties of rock.

To aid in interpretation of hydroseismograms, it is important that data on the strain wavefield also be collected at, or very near, the well. Broadband seismograms can be used to infer dynamic strain under certain assumptions about wavefield polarization and phase velocity; these assumptions often hold for long-period surface waves, but analyzing body waves can be problematic. A borehole strainmeter is an ideal instrument to compare with fluid pressure changes because it measures the strain wavefield directly, with high accuracy at both local and teleseismic distances. At the USGS well there is a Trillium 120 broadband seismometer at the surface (station GM.IWM01), and there is a borehole strainmeter nearby (station PB.AVN2); the strainmeter is located ~45 km away from the well, making it suitable for comparison with broadband seismometer records of teleseismic waves. The 2017 M_w8.2 Tehuantepec, Mexico, earthquake (us2000ahv0) revealed how closely linked ground velocities and strain can be (Figure 3). The close agreement between ground velocity and dynamic strain seen with the Tehuantepec event is confirmation that this is a highly coherent, planar wavefield; but, as Agnew and Wyatt (2014) showed, it is not safe to assume these conditions hold always. With colocated strain measurements, there is no need to rely on the assumption of planar wave propagation with a single phase velocity: strains are measured directly rather than inferred from acceleration. There is an additional benefit in that it is readily apparent when the assumption of a linear, undrained response breaks down (e.g., Barbour, 2015).

Figure  3.  Correspondence between fluid pressures, ground velocities, and tensor strains at the Arbuckle well for the 2017 M_w8.2 Tehuantepec earthquake. (a) Map of epicenter and well, similar to Figure 2. (b) Zoom-in of Osage County, Oklahoma, with the location of strainmeter AVN2, and mapped faults from the Oklahoma Geological Survey. (c) Simplified diagram of general layer thicknesses and depths of seismic and strain observations. (d) Timeseries from the well, seismometer, and strainmeter. The Love wave shows prominently in the T and E_rt records around 550 seconds. W_L: Relative fluid levels in the well, in mm. T, R, Z: transverse, radial, and vertical velocities in the direction of wave propagation, in mm/s. E_rt, E_rr, E_tt: Shear, radial, and transverse extension, in parts-per-billion (10⁻⁹). E_rr+E_tt is areal strain.

DownLoad: Full-Size Img PowerPoint

Shalev et al. (2016a) express Equation (1) in terms of areal and shear strain. In a radial-transverse coordinate system relative to the direction of wave propagation, the change in pore pressure from horizontal deformation is:

where $\kappa_u$ and $\nu_u$ are the undrained bulk modulus and Poisson's ratio, respectively, and $N$ is the shear-strain coupling factor. E_rr, E_tt, and E_rt are the tensor strain components; we assume near surface boundary conditions, where the vertical strain (E_zz) is related to the areal strain (E_rr+E_tt) by the expression $E_{zz} = -\nu_u(E_{rr} + E_{tt}) /$$(1 - \nu_u)$. While it would be preferable to estimate these coefficients with the full, 3D strain tensor, the type of strainmeter we use can measure only the horizontal strain tensor. Additional moduli may be needed if the system is transversely isotropic, for instance, but for linear poroelastic materials,

where $\mu$ is the elastic shear modulus. In the case of materials with an irreversible damage response (e.g., Shalev et al., 2016a), both $N$ and $A$ are complicated expressions involving additional parameters; however, we assume that the stress perturbations from our study events cause only reversible deformation.

To investigate whether shear coupling causes dilatant deformation in the Arbuckle reservoir, we supplemented the original teleseismic dataset (Barbour et al., 2019) with additional events that caused dynamic pressure changes, updated the magnitude-distance scaling relationship, and inverted the pressure and strain timeseries to solve for the volumetric and shear strain coefficients, $M$ and $N$. To prevent bias associated with the spatial separation of the strainmeter and the well, which introduces an apparent phase shift because of differences in arrival times, we use cross-correlation to align the strain timeseries with the ground velocity timeseries (see Methods). These are relatively small but detectable time shifts, so we restrict our analyses to events with both seismic and strain records (Table 1).

Table  1.  Teleseismic events used to estimate strain coupling ${\rm{coefficients}}^{\dagger}$.

Comcat ID	Origin time (UTC)	Mw	Lat. (°N)	Lon. (°E)	Depth (km)	$\theta$ (°)	$D$ (km)	$\Delta h$ (mm)
us2000ahv0	2017-09-08 04:49:19	8.2	15.02	−93.9	47.4	174	2418	277
us2000ar20	2017-09-19 18:14:38	7.1	18.55	−98.49	48	185	2023	17
us2000d3km	2018-02-16 23:39:39	7.2	16.39	−97.98	22	183	2258	15
us1000gez7	2018-08-21 21:31:47	7.3	10.78	−62.91	146.2	136	4431	21
us60003sc0	2019-05-26 07:41:15	8	−5.81	−75.27	122.6	156	5203	61
ci38457511	2019-07-06 03:19:53	7.1	35.77	−117.6	8	260	1893	19
us60007idc	2020-01-28 19:10:24	7.7	19.42	−78.76	14.9	142	2585	103
Notes: $\dagger$ These are events with records of not only high-frequency fluid pressure changes, but also broadband ground velocity at the well (GM.IWM01) and dynamic strains at strainmeter PB.AVN2. $\theta$: back azimuth from well to source (degrees clockwise from North); $D$: hypocentral distance to the well ($36.727^\circ$N, $-96.532^\circ$E); $\Delta h$: peak (absolute) change in fluid level.

 | Show Table

DownLoad: CSV

3.   Results: Source Scaling, Step Changes, and Dynamic Strain Coupling

In revising the magnitude-distance scaling of Arbuckle reservoir pore pressure changes from teleseismic waves, we also included data from seismically active regions dominated by natural seismicity. Specifically, we use data from Barbour (2015) for stations that show an undrained response to teleseismic waves, namely stations B084 and B088, in Southern California, and stations B076 and B078 in central California. We also use data from Ohno et al. (1997) for station EDY on the Izu Peninsula, Japan. The new relationship, for moment magnitude M_w and hypocentral distance $D$, is found by linear mixed-effects regression to be

with a residual standard error of 0.23 (base-10 log scale). The term $\delta_s$ is a random effect to account for site-dependent bias:

Figure 4 shows the relationship between the observed pressure change and the model prediction (Equation (5)). The close agreement between this magnitude-distance scaling relationship and one for coseismic water level changes observed along the San Andreas fault in central California, $0.89{{M}} - 1.6\log{D}$ (Roeloffs, 1998), suggests that coseismic and dynamic fluid pressure changes from distant earthquakes are largely borne out of the same physical process.

Figure  4.  Revised magnitude-distance scaling of peak dynamic pore pressure changes from teleseismic waves. Left: Observed pressure changes compared to those from Equation (5). The predicted values for the triggering events in van der Elst et al. (2013) are calculated at the distances to the Arbuckle well. Right: Data coverage with model prediction contours (in Pa) without adjustments for site-dependent effects (i.e., Equation (6)).

DownLoad: Full-Size Img PowerPoint

One complicating factor for interpreting a shear response is that there is a (quasi)static response to surface waves (Figure 5) similar to those first studied by Roeloffs (1998) along the San Andreas fault; this is surprising since the static strain from the rupture itself is immeasurably small at teleseismic distances, and ground motion intensities are much less than for local earthquakes where static strain is expected.

Figure  5.  Example of a static pressure response of the Arbuckle reservoir to passing surface waves from the 2017-09-08 M_w8.2 Tehuantepec earthquake (see also Figure 3). After the Rayleigh wave, average pressures inside the Arbuckle reservoir are seen to be higher than before the surface waves arrived.

DownLoad: Full-Size Img PowerPoint

There is considerable uncertainty as to the source of both static pressure increases and decreases caused by long period surface waves. The general understanding of this phenomenon is limited, but the most plausible explanations are related to changes in the local properties of the reservoir. For instance, Brodsky et al. (2003) proposed that static shifts in fluid level are related to change in the local permeability, and various models exist which explain the causes of both increases in decreases in pressure (Doan and Cornet, 2007; Manga et al., 2012; Shalev et al., 2016b); however, such a model does not necessarily permit both increases and decreases at the same well. Further, this physical explanation would most likely show a relationship to seismic energy density (e.g., Wang CY and Manga, 2010), but the static changes at the Arbuckle well show no such correlation (Figure 6). It is also challenging to understand how the local pore pressure regime could re-equilibrate so rapidly — tens to hundreds of seconds — and how this would be sustained over hundreds or perhaps thousands of events (geologic timescales). An alternative explanation is that the static changes are related to irreversible compaction or dilatancy effects (e.g., Shalev et al., 2016b). With respect to the dynamic response, cyclic deformation of preexisting cracks with low aspect ratios would also promote an azimuthally dependent response: such cracks, when aligned with the maximum principal stress, will tend to close first and promote an elastic and reversible anisotropy, assuming that contact stiffness increases with applied load. One future test of these models would be to monitor for changes in microseismicity rates within the reservoir over time, in the vicinity of pressure shifts.

Figure  6.  Dynamic and static pressure changes as a function of model pressure according Equation (5). The Tehuantepec observations are highlighted for reference (see also Figures 3 and 5). Compared to the dynamic response, the static response shows no correlation in either sign or amplitude with scaling based on source magnitude and distance.

DownLoad: Full-Size Img PowerPoint

For all events with dynamic strain records (Table 1), we inverted Equation (3) to solve for the coupling coefficients with additional terms to account for tidal trends and static shifts in the Arbuckle reservoir pressure response (see Methods). Among these events, we observe distinct classes of pressure response due predominantly to harmonic strains (Figure 7a) and more impulsive strains (Figure 7b); there are also large early-time signals that are not well-explained by dynamic strain (Figure 7c). The latter type of response, with more impulsive signals, suggests a response with azimuthal dependence.

Figure  7.  Observed teleseismic waves in Arbuckle reservoir fluid pressure changes, modeled by high-frequency dynamic strain measurements at PB.AVN2 (e.g., Figure 3). The map shows the azimuthal distribution of events with strain records, which are labeled, as well as the direction of the maximum horizontal stress (θ_SHmax, from Alt and Zoback, 2017). The inset map shows the location of the 2016 M_w5.8 Pawnee earthquake — the largest injection-induced earthquake on record. The records are grouped by apparent response type due predominantly to: (a) harmonic strains (e.g., Rayleigh waves), (b) impulsive strains (e.g., Love waves), and (c) large, early-time impulsive signals not well represented by dynamic strains. M is moment magnitude in this figure.

DownLoad: Full-Size Img PowerPoint

The strain coupling coefficients and static changes are given in Table 2. The areal strain coefficients ($M$) vary from 0.579 to 10.5 GPa, with standard errors up to 0.389 GPa; and the shear strain cofficients ($N$) vary from 0.189 to 7.99 GPa, with standard errors up to 0.591 GPa. The ratio of shear strain to areal strain coefficients ($N/M$) is a useful quantity to measure (see Methods) and here it shows less variation that the individual coefficients, ranging from 0.238 to 1.09. The static changes are robustly determined for all events, ranging from −33 to 159 Pa with standard errors up to ~9.5 Pa; static changes may be both positive and negative.

Table  2.  Areal and shear strain coupling coefficients (Equation (3)) for teleseismic events with dynamic strain records (Table 1).

ID	Step (Pa)	$M$ (GPa)	$N$ (GPa)	$N/M$	$\sigma_{\text{Step}}$ (Pa)	$\sigma_M$ (GPa)	$\sigma_N$ (GPa)	$p(M)$	$p(N)$
us2000ahv0	$159$	4.74	1.92	0.41	9.53	0.11	0.16	—	—
us2000ar20	$36.5$	4.14	1.82	0.44	1	0.18	0.25	—	0.01
us60007idc	$18.9$	1.03	1.02	0.99	4.88	0.11	0.11	0.1	0.52
ci38457511	$-16.7$	0.58	0.19	0.33	1.36	0.04	0.04	—	0.15
us60003sc0	$-25.6$	7.32	7.99	1.09	1.75	0.39	0.59	0.85	0.12
us1000gez7	$-30.9$	10.5	2.5	0.24	0.63	0.16	0.26	—	—
us2000d3km	$-33.1$	1.86	0.76	0.41	0.92	0.11	0.15	—	0.23
Notes: $\sigma$: standard error of the coefficient. $p(\cdot)$: p-value of F-tests on the coefficient; missing values represent probabilities below machine precision ($2 \times 10^{-16}$); the values for the step estimates are all below this threshold.

 | Show Table

DownLoad: CSV

To better constrain the average strain coupling behavior of the Arbuckle reservoir, we re-fit the pressure timeseries with a single inversion. First, all timeseries are fit with a weighted least squares regression procedure (see Methods), and compared with a linear mixed-effects regression. There is little difference between those results (denoted as ‘LM’ and ‘LME’ in Table 3, respectively). However, given that we observe poor agreement from the individual results for events ci38457511 and us60007idc (see Figure 7c), we further repeat the weighted least-squares regression after removing those problematic events from the dataset (result denoted as ‘LM (subset)’ in Table 3). This culling procedure is done only to confirm that the results are not overly sensitive to the few somewhat atypical observations.

Table  3.  Strain coupling coefficients for all events (Table 2).

${\rm{Fit}}^\dagger$	$M$ (GPa)	$N$ (GPa)	$N/M$	$\sigma_M$ (GPa)	$\sigma_N$ (GPa)
LME	3.56	1.28	0.36	0.04	0.05
LM	3.75	1.37	0.37	0.04	0.05
LM (subset)	4.74	1.94	0.41	0.04	0.06
Notes: $\dagger$: the fitting approach. LME is linear mixed effects regression, LM is ordinary least squares, and LM.sub is similar to LM but with poorly fitting events removed (see text, and Figure 7c).

 | Show Table

DownLoad: CSV

4.   Anisotropic Response of the Arbuckle Reservoir

We tested whether including shear strains in the fit is justified by computing an analysis of variance (ANOVA) table between a model that excludes shear strain — the null model — and a model that includes shear strain. We then evaluated the relative improvement (reduction) in variance with an F-test. For this comparison, the F statistic is 503.1, with a p-value less than machine precision ($2 \times 10^{-16}$). We conclude that it is appropriate to reject the null model in favor of the model that includes shear strain coupling. Taking the ‘LM.sub’ result as the fit most representative of the reservoir, the Arbuckle reservoir shows the following average strain-coupling behavior:

At this value of $N/M$, the value of Skempton's $A$ coefficient varies approximately linearly as a function of $\nu_u$ within the range for typical reservoir rocks ($\nu_u = [0.28,\ 0.33]$, Detournay and Cheng, 1993; see Methods). In this same range of $\nu_u$ the Arbuckle reservoir shows a response consistent with

This is notably smaller than 1/3, the theoretical value for isotropic materials; it is also similar to the value reported by Wang HF (1997) for limestone (0.22), and consistent with the lower end of the differential stress dependence reported by Lockner and Stanchits (2002) for sandstone (Figure 8). The poroelastic response of the Arbuckle reservoir includes a sensitivity to shear strain.

Figure  8.  Anisotropic poroelastic response of the Arbuckle wastewater disposal reservoir. Contours of the shear-to-areal strain coupling ratio (N/M) are shown as a function of the undrained Poisson’s ratio (ν_u) and Skempton’s coefficient A. The range of ν_u for typical reservoir materials is shown from Detournay and Cheng (1993). The range of A represents the stress dependence for Berea sandstone (Lockner and Stanchits, 2002) at its measured value of ν_u, 0.33. The dot is the value reported for Indiana limestone and tested on the observed, coseismic water-level response from a moderate earthquake on the San Andreas fault (Wang HF, 1997).

DownLoad: Full-Size Img PowerPoint

Kroll et al. (2017) found a Skempton's coefficient of $B = 0.58$ for the Arbuckle reservoir. This value, based on static pressure changes associated with the 2016 M_w5.8 Pawnee earthquake that were observed in similar Arbuckle monitoring wells, may apply to the Arbuckle reservoir at our monitoring well since the earthquake occurred ~50 km away (see inset map in Figure 7). Unfortunately, we do not have measurements of the Pawnee earthquake so we cannot confirm their preferred value; but, we note that shear strain coupling was neglected in their analyses which makes it possible that this value is biased high since, as Wang HF (1997) shows, shear coupling can have a pronounced effect on the spatial distribution of coseismic pore fluid pressure changes. Nonetheless we can use this estimate along with the laboratory value of $B \approx 0.45$ from Lockner and Stanchits (2002) to relate $N$ directly to the undrained shear modulus. Taking $A = 0.24$, Equation (4) can be expressed as:

which may then be related to other elastic properties known to vary spatially, such as S-wave velocity (e.g., $\mu = \rho V_S^2$). Based on Equations (8) and (D2), and the range of $\nu_u$ shown in Figure 8, we find:

These values are consistent with general poroelastic constants (Rice and Cleary, 1976; Detournay and Cheng, 1993; Wang HF, 2000).

Not only have these teleseismic data revealed a shear strain coupling effect in the poroelastic response of the Arbuckle reservoir, but they also show systematic, azimuthal variations that are most obvious in both the shear-to-areal ratio and the size of the static response (Figure 9). Ratios for individual events are clustered around the value from the full inversion (i.e., Equation (8)) except in the range of azimuths between those of shallow, mapped faults (Kolawole et al., 2019) and S_Hmin, the direction of the minimum horizontal stress (i.e., Alt and Zoback, 2017); in this region the ratios are at least twice as large. In contrast, while the absolute size of the static change is generally limited to < 40 Pa with no clear trend in sign (pressure increase or decrease), the 2017 M_w8.2 Tehuantepec earthquake arrived in the direction of S_Hmin and generated a step increase more than 4 times the average. Interestingly, the one event arriving in the direction of S_Hmax (event ci38457511) shows average behavior.

Figure  9.  Azimuthal variation in (a) the ratio of shear and areal strain coupling coefficients, and (b) and static pressure change. The vertical filled regions show the range of directions for faults mapped at the Arbuckle/basement contact (‘K19’: Kolawole et al., 2019) and for principal horizontal stress directions (Alt and Zoback, 2017, see also Figure 7). The blue ticks at the top show the approximate arrival direction from the 2010 M_w8.8 Maule and 2011 M_w9 Tohoku earthquakes, which van der Elst et al. (2013) report to have dynamically triggered earthquakes in Oklahoma. The 2017 M_w8.2 Tehuantepec event (largest symbol) arrived in the direction of the minimum horizontal stress, and caused both the largest dynamic and static response — a pressure increase.

DownLoad: Full-Size Img PowerPoint

Even though this dataset is relatively limited in coverage, the azimuthal dependence of the coupling coefficients suggest that the source of the shear strain coupling is not necessarily related to bulk properties. Laboratory observations of the poroelastic response of typical reservoir rock under realistic loading conditions show evidence for a differential stress dependency in anisotropy (Lockner and Stanchits, 2002; Lockner and Beeler, 2003; Wong, 2017), but it is difficult to to reconcile the sharp contrasts in values across tectonic indicators such as fault orientations and stress directions with such models. And, at $N/M$ ratios of ~1 or more (e.g., Figure 9a), the reservoir rock would have extreme anisotropy (e.g., Figure 8) that is likely to manifest in seismic velocities (Hudson, 1981; Cheng, 1993); this is not supported by recent studies that have documented relatively smooth variations in shear-wave splitting directions (Cochran et al., 2020) and crustal stress orientations (Alt and Zoback, 2017) in Oklahoma.

Instead, these data suggest an influential role for fluid-filled fractures. If the dynamic deformation during seismic waves modulates fracture apertures, and the pre-existing fracture pattern is relatively uniform, there would likely be an azimuthal variation in response, since the rock would be relatively compliant in the direction normal to fractures compared to the direction parallel to them. In other words, seismic waves might be effective only at modulating fracture apertures within a relatively narrow range of arrival directions, with that range depending on the state of stress and the predominant fracture pattern. Furthermore, the strong response seen in the direction of S_Hmin but not S_Hmax suggests that the static response is not due to fracture creation or growth, which, in this underpressured, strike-slip/normal faulting stress regime (S_v $\approx$ S_Hmax > S_Hmin), would occur roughly along a vertical plane striking parallel to S_Hmax. The dynamic stresses and pore pressure changes associated with these surface waves are relatively small, but might be large enough to mobilize particles that could either prop open or close off pre-existing fractures (e.g., Brodsky et al., 2003; Elkhoury et al., 2011). We also may not yet rule out the possibility of undrained consolidation (e.g., Wang CY and Manga, 2010) or dilatancy effects (e.g., Shalev et al., 2016b), although it is challenging to understand how these phenomena would show azimuthal variability. With the data at hand, we cannot determine whether the anisotropy is intrinsic (e.g., due to formation structure or bedding), stress induced, or both; this is the subject of future work.

5.   Conclusions

Considering the azimuthal, anisotropic response measured in the Arbuckle reservoir, and increasing evidence that basement rooted faults extend into the Arbuckle Group (e.g., Firkins et al., 2020), the primary implication for the triggering mechanisms of the events studied by van der Elst et al. (2013) is that a more complex interaction between dynamic stresses, pore pressure changes, and subsurface structure and rock properties is at play. Further, the details of static stress changes affecting earthquake-to-earthquake interactions should also consider the potential for coupling between pore pressure and shear deformation (e.g., Wang HF, 1997). These results also have implications for numerical simulations of fluid injection. With injection into a layered system like the Arbuckle reservoir, shear stress concentrations develop where there are strong contrasts in poroelastic properties (between adjoining layers). Additional differential stresses develop at seismogenic depths because of rapid reductions in crustal permeability with depth. As formation pressures in the Arbuckle continue rising due to wastewater disposal (Peterie et al., 2018; Barbour et al., 2019; Ansari et al., 2019), it may no longer be sufficient to neglect coupling between shear stress and pore fluid pressure changes, especially when considering how wastewater disposal causes failure stresses to develop at seismogenic depths.

Data and Resources

Seismic data from the monitoring station are accessible through Incorporated Research Institutions for Seismology Data Management Center (IRIS DMC) with station code IWM01, and network code GM. Strain data are from the Network of the Americas borehole strainmeter network, which is also accessible through the IRIS DMC (https://ds.iris.edu/ds/nodes/dmc/earthscope/pbo/) with station AVN2 and network code PB. Strain related material is based on services provided by the GAGE Facility, operated by UNAVCO, Inc., with support from the National Science Foundation (NSF) and the National Aeronautics and Space Administration under NSF Cooperative Agreement EAR-1724794. Fluid level data from the monitoring station are accessible in near-realtime through the station's NWIS site: https://waterdata.usgs.gov/nwis/uv/?site_no=364337096315401. Mapped faults in Oklahoma are from the Oklahoma Geological Survey report OF2-2016.

Websites were last accessed January 2021.

Acknowledgments

We thank David Lockner, Colin Pennington, the Editor Hongfeng Yang, and the three anonymous reviewers for providing thorough reviews that helped strengthen this work. We also thank Mai-Linh Doan for sharing a simple solution to the relationship between ground velocity and strain, and our colleagues at the U.S. Geological Survey Oklahoma Water Science Center for helping to maintain a high-quality Arbuckle reservoir observatory. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. government.

Appendix: Detailed Methods

A.   Estimating strain coupling

To invert for shear and areal strain coupling coefficients, we use weighted least-squares regression to fit the fluid pressure data to the areal strain, the Hilbert transform of the areal strain, the shear strain, and the Hilbert transform of the shear strain. The Hilbert-transformed series are used to stabilize the inversion because they are shifted relative to their original series by $\pi/2$ radians; this especially helps in fitting the approximately harmonic Rayleigh waves, for instance. We also fit a second order polynomial, and a Heaviside step function aligned at the time of the peak dynamic pressure change; these account for tidal trends and step changes in pressure, respectively. The strain coefficients are then found by summing the coefficients for each quantity in quadrature; their standard errors are computed with Taylor series expansion:

where $a$ is the strain regression coefficient (with standard error $\sigma_a$, $b$ is the Hilbert-transformed-strain coefficient (with standard error $\sigma_b$), and $c = \sqrt{a^2 + b^2}$. The modeled pressure changes shown in Figure 7 represent the best fitting coupling coefficients applied to the areal and shear strain series along with the step change estimate.

As an attempt to minimize the effects of overfitting in the case of low signal-to-noise ratios, the weighting scheme used in the regression ($W$) is set by the gamma distribution function:

where $\Gamma$ is the gamma function, $\gamma$ is the lower incomplete gamma function, $x$ is the quantile in consideration, $\alpha$ is the shape parameter, and $\beta$ is the scale parameter (see Abramowitz and Stegun, 1972, ch. 6). Here, $x = |\Delta p|$, the absolute value of the seismically induced changes in Arbuckle reservoir fluid pressure. We use $\beta = 1$, and $\alpha = n/l$, where $n$ is the number of observations for a particular event and $l$ is a constant that determines the downweighting cutoff. We use $l = 100$ to ensure a conservative cutoff. For instance, if the length of the series is 1000 samples, and the maximum pressure change is 200 Pa, this scheme reduces the influence of pressure variations below 20 Pa, with values of 10 Pa having half the weight as at 200 Pa. Decreasing $l$ for the same $n$ gives a higher pressure cutoff.

B.   Pre-processing of strain data

Prior to inversion, the linear gauge strains from the Gladwin-style four-component tensor strainmeter ($\boldsymbol G$, a $4 \times n$ matrix of $n$ strain observations from 4 strain gauges) are transformed to horizontal tensor strains ($\boldsymbol E$, a $3 \times n$ matrix) using the following coupling equation:

where ${\boldsymbol{C}}$ is determined by tidal calibration (e.g., Hodgkinson et al., 2013, see AVN2 station page):

The resulting matrix of tensor strains contains the areal strain and two engineering shear strains in an east-north coordinate system. With knowledge of the station-to-event azimuth $\theta$ (or back-azimuth), a coordinate transformation can be applied to rotate $E$ into a radial-transverse coordinate system ($\boldsymbol E'$) (e.g., Agnew and Wyatt, 2014):

C.   Relationship between ground velocity and strain

The areal strain for a Rayleigh wave in a homogenous halfspace can be calculated from the spatial gradients in the time-derivatives of the displacement functions (Ben-Menahem and Singh, 1981, e.g., their Eq. (3.105)) and it may be shown that this strain is proportional to the Hilbert transform ($\mathcal{H}$) of the vertical ground velocity (neglecting modulus-frequency dispersion):

where $V_P$, $V_S$, and $V_R$ are the P, S, and Rayleigh wave speeds, respectively. As discussed in the text, we align $\boldsymbol E'$ to the broadband ground velocity to minimize bias due to differences in traveltime between the well and the strainmeter. By way of the relationship in Equation (C1), we calculate the shift between the areal strain and the Hilbert-transformed vertical velocity, and apply the same shift to the other tensor strain components. Using least squares regression, we find an optimal value of $c = 1.89 \times 10^4 \;\text{m/s}$ for the aligned tensor strain and velocity timeseries. Figure 10 shows this best-fitting scaling coefficient applied to all records in this study with broadband seismic data, including the Alaska event shown in Figure 2.

Figure  10.  Broadband ground velocities scaled to strain using the relationship in Equation (C1) compared to Arbuckle reservoir fluid level timeseries. Plots are similar to the one shown in Figure 2. M is moment magnitude in this figure.

DownLoad: Full-Size Img PowerPoint

D.   Ratio of poroelastic strain coupling

Following Equation (3), the ratio of the shear and areal strain coupling coefficients is

and, from linear poroelasticity (Wang HF, 2000) the ratio of shear modulus to undrained bulk modulusis

Thus, Equstion (D1) can be simplified as: