A physics-informed optimization workflow to manage injection while constraining induced seismicity: The Oklahoma case

1 Introduction

In response to the rapid increase of the rate of seismicity in 2015 in central and northern Oklahoma, the regulators decided in 2016 to reduce the total volume of injected brine (waste-product of the hydrocarbon industry activities, notably shale gas production) to less than 40% of the 2014 total volume (Ellsworth, 2013; Walsh and Zoback, 2015). Late 2016, large magnitude earthquakes were recorded, and in April 2017 all the wells were shut in. In 2016–2017, the detailed causal relationship between injection operations and seismicity was unclear. In contrast, it is now well recognized and documented that the surge of the number of earthquakes in the central and northern Oklahoma can be attributed to the high-volume injection of brine in the subsurface (Keranen et al., 2014; Weingarten et al., 2015).

Multiple models have already been deployed to study either 1) the physical mechanisms at the origin of the induced events at Oklahoma (Norbeck and Horne, 2016; Goebel et al., 2017a; Goebel et al., 2017b; Johann et al., 2018; Dempsey and Riffault, 2019) or 2) to predict the return of the seismicity increase to a lower background rate after the stop of the injection activities in 2017 (Langenbruch and Zoback, 2016; Goebel et al., 2017a; Goebel et al., 2017b; Langenbruch et al., 2018; Zhai et al., 2019). Up to now, these modelling strategies always involved simplifications either at the level of the simulations of the flow throughout the porous media (Goebel et al., 2017a; Goebel et al., 2017b; Zhai et al., 2019) or at the level of the computations of the induced stress development and seismicity changes (Langenbruch and Zoback, 2016; Norbeck and Horne, 2016; Johann et al., 2018; Langenbruch et al., 2018; Dempsey and Riffault, 2019). In order to gain confidence in the predictive power of any modelling strategy, a thorough model validation against observed data is required. This step has largely been overlooked in existing modelling approaches of induced seismicity at Oklahoma which generally rely on simplified sensitivity analysis.

In this study, a novel modelling approach is outlined that is tailored to honor as much as possible all available pre-existing knowledge at injection sites in the central Oklahoma area (for simplicity the area is simply referred to as “Oklahoma” in the remainder of the manuscript, see Figure 1). Available knowledge consists of data on local geology, flow, mechanical and seismicity response. The model parameters of our seismological model are conditioned with the observed seismicity in order to best fit the full spatio-temporal history of induced earthquakes. It is only after completing this data assimilation procedure that injection strategies can be optimized to maximize the volume of injected brine while minimizing induced seismicity. For optimizing injection strategies, we deployed a constrained optimization workflow that maximizes injection volumes under a fixed limit of allowed seismicity rate. Optimization results show with a different spatio-temporal injection history, the rapid growth of the seismicity rate in 2015 can be successfully prevented while the total volume of injected brine can be maintained. It indicates that injection operations and mitigation measures for induced seismicity greatly benefit from optimization of spatio-temporal injection strategies as seismic risks can be reduced under continuing injection operations.

FIGURE 1

FIGURE 1. Maps of the Oklahoma state in United States showing: (i) the area of interest (dashed contour in top-right map), (ii) the locations of the Arbuckle wastewater disposal wells (black dots in bottom-centre map), (iii) and the M_w3+ earthquakes which nucleated in the basement between 3 and 7 km depth and from 1995 to 2017 (red stars in the bottom-centre map).

2 Methodology

2.1 Flow simulations

The two main identified mechanisms that cause induced seismicity at Oklahoma are: 1) the slow diffusion of elevated pore pressure from the high-permeability Arbuckle aquifer (locus of the brine injections) to its low-permeability basement at the nucleation depth of the induced events (e.g., Langenbruch and Zoback, 2016; Norbeck and Horne, 2016; Dempsey and Riffault, 2019); 2) the change in total stress at this depth in the basement that is caused by poro-elastic loading due to the inflation of the rock volume associated with the increase in pore pressure (e.g., Goebel et al., 2017a; Goebel et al., 2017b; Zhai et al., 2019).

The first modelling effort consists of assessing the historical spatio-temporal distribution of pore pressure changes caused by the high-volume injection of brine in the subsurface at Oklahoma. As in each step of our modelling strategy, this step aims to incorporate all the pre-existing available knowledge on local geology, flow, mechanical and seismicity response at Oklahoma (Johnson, 2008; Faith et al., 2010; Holland, 2013; Keranen et al., 2013; Hincks et al., 2018; Pei et al., 2018). The computation of the pore pressure changes (Figure 2) is performed using the open source OPM-flow simulator (Open Porous Media Initiative, 2020; version 2020.04). OPM-flow is a three-dimensional three-phase fully-implicit reservoir simulator with realistic representation of reservoir geological properties. Wells are incorporated with the standard models used in the commercial simulators from the oil and gas industry. Multiple options are available, as black-oil, thermal, aquifer, CO₂; in our specific case we modelled brine. The reader is referred to Rasmussen et al. (2021) for a complete description of OPM-flow.

FIGURE 2

FIGURE 2. Flow simulation—visualization via ResInsight [open source visualization software, ResInsight (Computer Software), version 2020.10, part of Open Porous Media Initiative, Ceentron Solutions, 2020] of the pressure field (changes in pore pressure since the start of injection; in bars) at the end of 2015. Vertical exaggeration x5.

The flow model is populated by properties based on current understanding of the geological and hydrogeological settings at Oklahoma. Three layered-box models have been deployed, each of them including the same thickness for the Arbuckle aquifer and its basement but with different permeabilities and porosities (Table 1). Only the results of the flow simulation with a permeability of 50 and 0.05 mD for the Arbuckle aquifer and its basement, respectively, are discussed in this report (i.e., Model 1 in Table 1). This permeability contrast yields the most likely changes in pressures in the Arbuckle and its basement and is preferred to best match observations.

TABLE 1

TABLE 1. Geometry and model parameters of the flow simulation models.

The OPM-FLOW simulation model honors the geographical location and depths of all the wells according to the Oklahoma Corporation Commission. It includes 220 wells for central Oklahoma, the focus area for the present study, and for each of the well injectors the historical monthly injection rate from January 1995 to January 2018 is used as input (Figure 3). Note here that central and western Oklahoma can be considered as isolated and independent compartments in terms of flow (see Zhai et al., 2019); the same modelling complexities are attached to both areas of interest and solely focusing on central Oklahoma is considered sufficient to demonstrate the capabilities of the newly developed modelling framework. The injection depth at the Arbuckle level and distance to the basement is included in our modelling strategy and has been reported as an important parameter correlated with the occurrence of the induced events (Hincks et al., 2018). The run-time of one full OPM-FLOW simulation from January 1995 to January 2018 takes less than 20 min, which is rather fast for such a complex flow simulation. The low computation time makes the simulations suitable for an ensemble-based approach that is used for the optimization. This optimization is the last step of our modelling approach (see Section 2.6).

FIGURE 3

FIGURE 3. Historical yearly field-wide injection rate.

2.2 Induced stress in the basement

The objective of this modelling step is to assess the spatio-temporal development of induced stresses at the earthquake nucleation depth. Most of the events of the catalog compiled by the Oklahoma Geological Survey are located at a depth between 3 and 7 km (McGarr and Barbour, 2017; Pei et al., 2018). Considering an average depth uncertainty of 2 km, an average earthquake nucleation depth of 4 km depth is considered for our analysis. The two main mechanisms that were identified to control induced stress development and earthquake nucleation at Oklahoma are included: 1) the decrease of the effective normal stress along the basement faults by pore pressure diffusion and 2) the changes in total stresses induced by poro-elastic effects associated with pressure increase and volumetric expansion of rock.

Our modelling approach assumes matrix-dominated flow in the Arbuckle aquifer and fracture-dominated flow in the basement. More specifically highly permeable faults of the basement are assumed to hydraulically connect the base of the porous Arbuckle aquifer to the earthquake nucleation depth. Simple analytical calculations including typical fault properties (diffusivity, porosity, thickness) show that the diffusion time throughout the highly permeable faults-channels is on the order of at most a few days (Zhai et al., 2019). Therefore, one can further assume that the changes in pore pressure modelled by OPM-FLOW at the base of the porous aquifer are representative of the changes in pore pressure at 4 km depth experienced by the basement faults. These changes in pore pressures relatively to the start of injection and along a horizontal plane at 4 km depth are shown in Figure 4. The increase in pore pressure is at first (i.e., 2010) localized at the center of the area of interest and then starts to progressively migrate towards the North-West.

FIGURE 4

FIGURE 4. Changes in pore pressure (in MPa) at 4 km depth along the basement faults.

From the spatio-temporal changes in pore pressures at 4 km depth (Figure 4), the changes in the normal effective stress can be directly determined. Analogous to the analysis by Zhai et al. (2019), vertical strike-slip faults with a unique preferred fault strike azimuth of 50° are assumed to be ubiquitous in the basement. This faulting regime and the fault orientation is consistent with the statistical analysis of earthquake focal mechanisms, in situ stress analysis, field mapping, and 3D seismic interpretation at Oklahoma (Holland, 2013; Alt and Zoback, 2016; Kolawole et al., 2019; Qin et al., 2019; Firkins et al., 2020). Principal stresses are assumed to be spatially uniform around Oklahoma with a maximum horizontal stress azimuth oriented at 85°.

Changes in total stress imposed by the volumetric changes of the rock volume is determined using the in-house mechanical simulator MACRIS (Mechanical Analysis of Complex Reservoirs for Induced Seismicity, van Wees et al., 2019; Candela et al., 2019). The main advantage of MACRIS is that it is a mesh-free simulator, i.e., it does not need construction of a dedicated grid for the geomechanical analysis. MACRIS directly takes the grid of the flow simulation as input. In the present case, the grid consists of the 3D pressure fields computed by OPM-FLOW at a yearly sampling rate. Each grid block of the flow simulation is considered as an inflating/compacting nucleus of strain (or center of inflation/compaction, see Mindlin 1936; Geertsma, 1973; Okada, 1992). The contribution of each of these nuclei is integrated to compute the poro-elastic stress change at 4 km depth in the basement. The Barnes-Hut algorithm (Barnes and Hut, 1986) is used for re-discretizing the initial flow grid. The purpose of using this algorithm is twofold: 1) to cluster the nuclei of strain close to the 4 km depth horizontal plane in order to increase the spatial stress resolution, and 2) to shorten the computation time. The approach has been validated by comparison with relatively slow finite-element numerical computations in a previous study (van Wees et al., 2018). The poro-elastic normal and shear stress changes acting on the faults at 4 km depth can thus be calculated using MACRIS. Coulomb stress changes can be evaluated by considering shear stresses and effective normal stresses:

$S=\tau -\left[\mu -\alpha \right]{\sigma }_n^{\prime }(1)$

where $\tau$ is the shear stress acting along the fault plane, ${\sigma }_n^{\prime }$ is the effective normal stress (which include the direct effect of the pore pressure increase intra-fault), $\mu$ is the coefficient of fault friction and $\alpha$ is a constitutive parameter (zero in this study). The Coulomb stress rate $\dot{S}$ along the horizontal plane at 4 km depth is presented in Figure 5. High Coulomb stress rates are first (i.e., 2010) localized in the center of the area of interest (hereafter called Central area), and then progressively migrate towards the North-West (hereafter called North area). In 2015, Coulomb stress rates are highly localized in the North area, and following the new measures imposed by the regulator, the Coulomb stressing rate starts to decrease and delocalize in 2016.

FIGURE 5

FIGURE 5. Coulomb stress rate fields (in MPa/Year) at 4 km depth along the basement faults. The magenta and green dashed circles in the top-left figure indicates the North area and Central area respectively.

2.3 Induced seismicity

The traditional Coulomb failure model predicts that whenever the Coulomb stress reaches a threshold value, an earthquake is generated (e.g., Ader et al., 2014). Assuming a population of faults on which the pre-stresses are uniformly distributed up to the threshold value, the Coulomb failure model depicts direct proportionality between the seismicity rate and the Coulomb stress rate. For example, during any arbitrary stressing history, the Coulomb failure model predicts an instantaneous reduction or rise of seismic events as soon as the Coulomb stress starts to decrease/increase. This prediction is not in agreement with the observed seismicity in Oklahoma since only a few observed events have been recorded in the year 2009 (cf. Section 2.4). The Coulomb failure model would predict a high seismicity rate that is linearly related to the high Coulomb stress rate in 2010 (Figure 5).

One shortcoming of the Coulomb failure model is that it does not honor the frictional constitutive behavior of faults. Laboratory data show that the timing of dynamic instability of faults depends on initial stress conditions, fault properties and applied stress (Dieterich and Kilgore, 1996). Rate-and-state friction laws have been established in order to reproduce these laboratory observations (see Marone, 1998 for a review). More specifically, the rate-and-state friction laws reproduce the fact that the onset of frictional sliding is a non-instantaneous time-dependent process (as opposed to the instantaneity assumption of the Coulomb model), which introduces a time-dependent failure mechanism for the generation of earthquakes. Now assuming a population of faults following a rate-and-state frictional behavior, and where the time-to-failure of the nucleation spots along the faults is uniformly distributed, Dieterich (1994) derived the following seismicity rate model:

$R_D=\frac{r_0}{\gamma {\dot{S}}_0}\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\frac{d\gamma }{dt}=\frac{1}{A{\sigma }_n^{\prime }}\left[1-\gamma \frac{dS}{dt}\right](2)$

and where $R_D$ is the seismicity rate, $\gamma$ is a state variable, $S$ is the modified Coulomb stress function defined in Eq. 1. The constant $r_0$ is the steady-state background seismicity rate at the reference stressing rate ${\dot{S}}_0$. $A$ is a dimensionless fault constitutive parameter.

Segall and Lu (2015) reformulated this seismicity rate equation to eliminate the state variable $\gamma$. They defined a normalized seismicity rate, relative to the background rate, as:

$R=\frac{R_D}{r_0}(3)$

The differential equation for $R$, derived from Eqs 2, 3, is:

$\frac{dR}{dt}=\frac{R}{t_a}\left[\frac{\dot{S}}{{\dot{S}}_0}-R\right](4)$

where $t_a=A{\sigma }_n^{\prime }/{\dot{S}}_0$ is the characteristic time delay for the earthquake nucleation process. This delay also corresponds to the time scale of the decay in aftershock rate from the main shock back to the background rate.

2.4 Declustering of the observed seismicity

The earthquake catalogue used for our study has been compiled by the Oklahoma Geological Survey. Our analysis is carried on with earthquakes larger or equal to M_w = 3; this minimum magnitude is above the completeness magnitude of the catalogue, and has been selected for sake of comparison with previous studies (e.g., Langenbruch and Zoback, 2016; Zhai et al., 2019). In addition, based on an average depth uncertainty of 2 km and to focus our analysis on the basement only the events nucleating between 3 and 7 km depth are selected. This depth range includes most of the events of the catalog. When this magnitude and depth filtering is applied to the raw earthquake catalogue, the maximum yearly rate of events at Oklahoma reaches roughly 400 events/year in 2015 (Figure 6).

FIGURE 6

FIGURE 6. Declustering effect providing a separation between mainshocks and aftershocks induced by injection activities Left: all events of the catalogue (grey line) and independent mainshocks (black line). Right: Bimodal distribution of time and space components (${\mathit{log}}_{10}T,{\mathit{log}}_{10}\mathcal{R}$) of the nearest-neighbor distance $\eta$. Each gray circle in these plots corresponds to an event in the seismicity catalogue. The location of the event in the (${\mathit{log}}_{10}T,{\mathit{log}}_{10}\mathcal{R}$) plane provides information about the time and space distance to the event’s parent. The threshold value ${{\mathit{log}}_{10}\eta }_0$ to separate the two modes is estimated used a 1D Gaussian mixture model applied to the logarithmic nearest-neighbor distances ${\mathit{log}}_{10}{\eta }_{ij}$ (Hicks, 2011). The red dashed diagonal lines depict the ${{\mathit{log}}_{10}\eta }_0$. On the right of the red dashed diagonal line, independent mainshocks are indicated, while events on the left of the line are clustered aftershocks.

One of the assumptions in Dieterich’s seismicity rate theory (Dieterich, 1994) that is subject of debate, is the lack of interactions between seismic sources. More specifically, Dieterich’s seismicity rate theory assumes aftershocks are directly triggered by the stress changes induced by the mainshock. Effect of stress interactions between aftershocks is not accounted for. To circumvent this potential shortcoming of Dieterich’s theory, we apply the declustering algorithm of Zaliapin and co-workers (Zaliapin et al., 2008; Zaliapin and Ben-Zion 2013; Zaliapin and Ben-Zion, 2016). This algorithm resolves complete triggering chains based on the relative space-time-magnitude distances, originally introduced by Baiesi and Paczuski (2004). Only the main ingredients of the approach are given below, but the reader is referred to the papers of Zaliapin and co-workers for more details.

The declustering algorithm of Zaliapin et al. (2008) is based on space-time-magnitude distances between earthquakes $i$ and $j$ defined as (Baiesi and Paczuski 2004):

${\eta }_{ij}=\{\begin{array}{c}{t}_{ij}{\left(r_{ij}\right)}^d{10}^{-b{m}_i},t_{ij}>0;\\ \infty ,t_{ij}\le 0.\end{array}(5)$

$t_{ij}$ is the event intercurrence time in years, which is positive if event $j$ occurred after event $i$ ($t_{ij}=t_j-t_i$); $r_{ij}\ge 0$ is the spatial distance between the earthquake hypocenters in kilometers; $m_i$ is the magnitude of event $i$; $d$ is the (possibly fractal) dimension of the hypocenters; and $b$ corresponds to the $b$-value of the Gutenberg-Richter frequency-magnitude distribution.

The nearest-neighbor distance for a given event $j$ is the minimum distance among ${\eta }_{ij}$ where $i$ goes over all earlier events in the catalogue. The event $i$ that corresponds to the nearest-neighbor distance is called the nearest-neighbor, or parent, of event $j$.

Zaliapin et al. (2008) proposed to consider the scalar distance $\eta$ in terms of its space and time components normalized by the magnitude of the parent event $i$ as:

$T_{ij}=t_{ij}{10}^{-qb{m}_i};{\mathcal{R}}_{ij}={\left(r_{ij}\right)}^d{10}^{-pb{m}_i};p+q=1(6)$

And now:

${\mathit{log}}_{10}{\eta }_{ij}={\mathit{log}}_{10}{T}_{ij}+{\mathit{log}}_{10}{\mathcal{R}}_{ij}.(7)$

For our analysis, we fix $b=1$, $d=1.6$, and $p=0.5$ following Zaliapin and Ben-Zion (2016). Zaliapin and Ben-Zion (2013) have demonstrated that the estimated cluster structure is stable with respect to the parameter values. Accordingly, our conclusion are expected to be non-sensitive to the precise parameter values.

Zaliapin et al. (2008) and Zaliapin and Ben-Zion (2016) have shown that observed seismicity generally presents a bimodal joint distribution of (${\mathit{log}}_{10}T,{\mathit{log}}_{10}\mathcal{R}$). In the present Oklahoma injection case, one mode corresponds to the independent mainshocks, whereas the other consists of clustered aftershock events located considerably closer in time and space to their parents. Figure 6 shows the 2D joint distributions of the rescaled time and space component ($T$, $\mathcal{R}$) of the nearest-neighbor earthquake distance $\eta$. For Oklahoma, after applying the declustering step, the maximum yearly rate of mainshocks reaches roughly 120 event/year in 2015 (see Figure 6).

2.5 Data assimilation

The objective of the data assimilation procedure is to determine the set(s) of model parameters [$r_0,{\dot{S}}_0$, ${\sigma }_s$] (see previous section for definition of symbols) in Dieterich’s (1994) seismicity rate theory (cf. Eq. 4) that give the best agreement between the observed seismicity rate and the computed rate.

The log-likelihood is defined as the logarithm of the probability that one specific model (with one specific set of model parameters) has generated the observed earthquake catalog. For each set of model parameters the posterior log-likelihood is calculated in order to rank the models (Ogata, 1998; Zhuang et al., 2012). The model with the highest log-likelihood is most likely to have generated the observed seismicity catalogue. For our non-homogeneous stationary Poisson process, and for a given time interval [t₀, t₁] and spatial area [x₀,x₁] × [y₀,y₁], the log-likelihood with respect to N observed earthquakes which have occurred at times t_i and locations (x_i,y_i) is defined by the following function:

$ll={\displaystyle \sum }_{i=1}^N\log \left(R_D\left(xi,yi,ti\right)\right)-\underset{t0,x0,y0}{\overset{t1,x1,y1}{\iiint }}{R}_D\left(x,y,t\right)dxdydt(8)$

Following a Bayesian approach, and because we do not have a prior preference for the shape of the distribution of each model parameter, we attribute a bounded uniform probability distribution for each model parameter. The Markov Chain Monte Carlo (MCMC) algorithm is used to condition these uniform distributions with the observed seismicity.

The start of the time period where the data assimilation procedure is deployed is 1st January 2009, start of the observed seismicity, ending on 31 December 2017. Note that we apply our analysis to seismicity rate only. Incorporation of seismic magnitude requires further steps in the approach. Ideally, a joint log-likelihood is applied to assimilate both seismicity rate and seismic magnitude. However, in order to properly model the magnitude of induced events, a much more complex modelling strategy is required combining the nucleation rate from Dieterich’s (1994) theory with, for example, the state of stress and energy available around each nucleation point (e.g., Noda et al., 2009; Schmitt et al., 2015; Dempsey and Suckale, 2016). As performed in previous studies at Oklahoma (e.g., Zhai et al., 2019), the nucleation rate could be combined with an arbitrary frequency-magnitude distribution that is uniform and constant in space and time. However, applying such a procedure would bring additional uncertainties and, for the current purpose of the study, little added value. Therefore, we focus the analysis on nucleation rate of seismic events alone.

2.6 Constrained optimization

After applying the data assimilation procedure, the best set of model parameters can be selected with which our predictive model is more likely to explain the observed induced seismicity rate. Using this calibrated predictive model and set of model parameters, the optimum injection strategy is determined. The aim is to prevent the peak of seismicity kicking off in 2014, while at least the same volume of injected brine is injected as has been historically reported for the area. The objective is to maximize the total cumulative field-wide injected brine volume. We include a threshold value for the seismicity rate as a constraint that cannot be exceeded. The objective (maximizing injected volume) and the constraint (minimizing seismicity rate) are expected to be conflicting. The identification of an injection strategy that meets both our objective and constraint is a complex task. The complexity is increased by the non-regular distribution of wells, and the time-dependency of the geomechanical-seismological response (linking pore pressure diffusion and earthquake nucleation) to changes in injection rate. We therefore adopt a numerical optimization approach that aims to solve the following formalized problem:

$\widehat{u}=\underset{u}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}\mathrm{x}}J\left(u\right)\mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}R\left(u\right)<R_{\mathit{max}};J\left(u\right)>J_{\mathit{min}}(9)$

The control vector $u$ contains the well injection rates which are constant during discretized time intervals, i.e., $u=\left[u_1^1,u_2^1,\dots ,u_{N_t}^1,u_1^2,\dots ,u_{N_t}^{N_w}\right]$ where the subscripts indicate the time interval, and the superscripts indicate the well. The total number of controls is expected to be very large in our scenario since there are more than 200 wells and multiple time intervals. $J\left(u\right)$ is the total volume of injected brine, and $R\left(u\right)$ is the maximum field-wide yearly event rate. $J_{\mathit{min}}$ and $R_{\mathit{max}}$ are the historical total cumulative field-wide injected brine volume from 1st January 1995 to 31 December 2017, and the threshold value of the field-wide yearly event rate during the same time period, respectively.

We use the Stochastic Simplex Approximate Gradient (StoSAG) estimation method where the gradients are used in an advanced method for constrained optimization (Chen et al., 2009; Fonseca et al., 2014, 2017). An ensemble of randomly chosen control vectors (well injection rates) is generated and the stochastic gradient, total volume of brine injected and yearly event rate are computed. Based on the stochastic gradient direction, the controls are updated, and a new ensemble of perturbed controls is regenerated. This iterative process continues until there are no more significant changes in the total amount of brine injected or controls, or when a maximum number of iterations is reached (Figure 7). As an initial estimate (i.e., first guess) for the controls we prescribe constant and equal rates for all wells, so that the cumulative injected volume is slightly lower than the historical total injected volume.

FIGURE 7

FIGURE 7. Constrained optimization scheme.

3 Results

3.1 Posterior model parameters of the unsteady Oklahoma fault system

The joint and marginal posterior probability distributions of the model parameters obtained after assimilation of the declustered catalogue over the period from 1st January 2009 to 31 December 2017 are well constrained (Figure 8). The median of the marginal posterior distributions of each model parameter define the best posterior estimates ($A=0.00027,r_0=18event/year,{\dot{S}}_0=$ 0.0004 MPa/year). With our modelling approach, the average effective normal stress is 44 MPa at 4 km depth, which leads to the best posterior estimate of the characteristic relaxation time of seismicity, $t_a=\frac{A{\sigma }_n^{\prime }}{{\dot{S}}_0}=29.7years$_.

FIGURE 8

FIGURE 8. Posterior probability distributions of the seismicity rate model parameters obtained after the data assimilation procedure with the MCMC search. Histogram-type plot: marginal probability distribution of each individual parameter $A$ (left), ${\dot{S}}_0$ (centre), $r_0$ (right). Cloud of points: joint probability distribution of each pair of model parameters $A$-${\dot{S}}_0$ (top-left), $A$-$r_0$ (bottom-left), $r_0$-${\dot{S}}_0$ (bottom-right). Note that the prior probability distribution (not shown on this figure) of each model parameter was bounded and uniform such as: $A$ (-): U($1e^{-5}$, $1.0$), ${\dot{S}}_0$ (MPa/year): U($1e^{-7}$, $1.0$), $r_0$ (event/year): U($0.001$, $2e^5$)—with U(a,b) is a uniform distribution between a and b.

The assimilation of observed events starts in 2009 when a significant earthquake activity started to be recorded. Consequently, the Dieterich (1994) seismicity rate Eq. 4 is solved assuming initial condition at steady state, that is $R\left(0\right)=$ 1 in 2009. However, the start of the human-induced perturbation of the Oklahoma fault system by high-volume brine injections goes back to 1995, and hence it is most likely that the background activity was not at steady state in 2009. In fact, defining $R^{\prime }=\beta R$ and by operating a change of variable, we can show that $R^{\prime }$ satisfies an equation similar to Eq. 4 and reads:

$\frac{d{R}^{\prime }}{dt}=\frac{R^{\prime }}{A{\sigma }_n^{\prime }/{\dot{S}}_0^{\prime }}\left[\frac{\dot{S}}{{\dot{S}}_0^{\prime }}-R^{\prime }\right](10)$

with the unsteady initial condition $R^{\prime }\left(0\right)=\beta$ in 2009, and the background stressing rate ${\dot{S}}_0^{\prime }={\dot{S}}_0/\beta$. Following Eq. 3, this solution corresponds to a background seismicity rate $r_0^{\prime }=r_0/\beta$. This mathematical derivation implies that the background stressing rate ${\dot{S}}_0^{\prime }$ and background seismicity rate $r_0^{\prime }$ were $1/\beta$ times lower in 1995, at the start of injection when the Oklahoma fault system is assumed to be steady, compared to the values of ${\dot{S}}_0$ and $r_0$ inferred for 2009. The fault constitutive parameter $A$ is expected to remain constant in time. However, the background stressing rate and seismicity rate inferred for 2009 should not be interpreted as real steady state background values (before the start of human-induced perturbations by brine injections), but should be understood as “apparent background values,” which are actually the reference values at the initial time of the analysis, here 2009.

3.2 Temporal pattern

A posterior ensemble of yearly event rate predictions can be visualized by randomly picking members from the posterior density distributions obtained with the MCMC search. This ensemble can be compared with the data. However, it should be noted that the modelling strategy does not yet include the intrinsic Poisson variability of earthquake occurrence. The observed declustered catalogue is considered here as one unique realization of a stochastic, non-stationary Poisson process. Each posterior member of the modelling approach is one particular model of the time-dependent seismicity rate underlying the Poisson process. For an appropriate assessment of the predictive performance of our models, the stochastic Poisson variabilities need to be accounted for. For each posterior member, which can be considered as the mean of an underlying Poisson distribution, multiple synthetic catalogues are generated where the likelihoods of the event location and timing are proportional to the event density (Zhuang and Touati, 2015). Figure 9 shows that the posterior predictions capture very well the temporal variation in observed seismicity. More specifically, our modelling strategy can predict:

(1) the relative seismicity quiescence from 2009 to 2013,

(2) the abrupt ramp-up of the seismicity rate starting in 2014,

(3) the fast decrease of the seismicity rate starting in 2016 with the new measures imposed by the regulators.

FIGURE 9

FIGURE 9. Comparison of the predicted seismicity histories with the data (black line). The model corresponds to 300 realizations (cyan) randomly drawn from the posterior density distribution obtained with the MCMC search during calibration. The mean of the models is indicated by the blue line. For each posterior member, 30 synthetic catalogues are generated in order to account for stochastic Poisson variabilities. The grey region indicates 95% of the distribution of the 9,000 synthetic catalogues when stochastic Poisson variabilities are accounted for. Top: cumulated number of events. Bottom: yearly event rate.

3.3 Spatial pattern

Predictions of spatial distribution of seismicity rate based on our modelling strategy can also be evaluated. Figure 10 show a comparison of the spatio-temporal distribution between observations and model predictions (i.e., event densities). The overall spatio-temporal evolution of the observed seismicity is well reproduced by our modelling strategy. More specifically, the observed and modelled seismicity both start to increase at the center of the area of interest (i.e., Central area) and then progressively migrate outwards to finally be concentrated in the North area.

FIGURE 10

FIGURE 10. Mean posterior fields of event density (/cell) at 4 km depth along the basement faults. The black dots indicate the observed events.

3.4 Towards an optimum injection strategy

Now that parameters have been adjusted and that our model best explains the historical observations, this calibrated model can be used for the optimization exercises. Multiple optimization experiments have been performed but in the present manuscript we focus on the description of one key experiment. For this specific key experiment, we seek to find the optimum injection strategy for each well for the last 4 years [2013–2017] leading to a total volume of injected brine that is at minimum equal to the total historic injected volume while keeping the yearly event rate below a certain threshold. This threshold is set to 70 event/year in order to avoid the peak of seismicity rate starting in 2014. More specifically, the injection rate can be updated each year of the last 4 years independently for each well, which corresponds to 880 controls. For each iteration, the number of perturbations is set to 100 to appropriately compute the stochastic gradients. The optimizer successfully converged towards an optimum solution which satisfies both constraints: 1) the total volume injected is quasi-identical to the historical total volume injected, and 2) the yearly event rate is equal to the threshold of 70 event/year (Figure 11).

FIGURE 11

FIGURE 11. Comparison of the cumulative field-wide injected volumes (up) and seismicity rates (down) for the 4 years of the optimization period. The horizontal dashed line (bottom graph) indicates the imposed threshold of 70 event/year. Note that the historical seismicity rate (blue curve in bottom graph) is similar to the mean posterior (blue curve) in Figure 9.

Instead of preferentially using specific wells for most of the injected volume as historically, the optimum scenario suggests a more spatially uniform distribution of the injected volumes (see Figure 12).

FIGURE 12

FIGURE 12. Comparison of the total injected volumes (m³) per wells. Note the range difference of each color scale in the bottom maps.

It is more interesting to identify any differences in the spatio-temporal pattern between the first guess scenario and optimum scenario. Indeed, the first guess scenario consists to the exact spatially uniform scenario where we prescribe constant and equal rates for all wells (see Figure 12). Therefore, any deviations from this first guess scenario should illuminate a non-uniform spatial pattern for the optimum scenario with preferential use of specific wells. It is difficult to draw conclusions solely based on the visual comparison of the event densities for the first guess and optimum strategies (see Figure 13). The spatial pattern seems identical for both strategies, except the overall higher event densities for the optimum scenario. However, and interestingly, plotting now the differences between the optimum and first guess event densities (see Figure 14), it appears clear that the optimizer preferentially used the wells of the North area. This optimum strategy can be explained by the fact that before the optimization period (that is before 2014) most of the wastewater have been injected in the Central area (see Figure 12) resulting in a very localized high event density in this area (see Figure 15). Therefore, the use of wells at the Central area has been penalized by the optimizer since additional stress changes in this area would have led to a yearly event rate which would have probably exceeded the threshold of 70 event/year.

FIGURE 13

FIGURE 13. First guess and optimum event densities (/cell) at 4 km depth along the basement faults. The black triangles indicate the locations of the injector wells.

FIGURE 14

FIGURE 14. Difference between the event densities of the optimum and first guess scenarios. The grid cells underlined in black correspond to the ones with a high cumulative event density before the optimization period (see Figure 15).

FIGURE 15

FIGURE 15. Cumulative event density (/cell) right before the start of the optimization period (end of 2013). Grid cells of the Central area with a number of events higher than 1.0 are underlined in black.

Figure 16 illustrates this last conclusion. Indeed, Figure 16 depicts, for the optimum scenario and one hypothetical scenario, typical stressing rate histories and relative seismicity rate histories for the Central area and North area. The hypothetical scenario illustrates the case where the injections and Coulomb stress rates would have been maintained in the Central area during the optimization period. In this hypothetical scenario the relative seismicity rate (already in a steady-state regime because of the high-injected volumes before 2014) would have been much higher than the optimum scenario, this last consisting in promoting injection in the North area.

FIGURE 16

FIGURE 16. Synthetic examples of the Coulomb stress and relative event rate evolutions for the North area and Central area (see Figure 5). The Central area (optimum) curve and North area (optimum) curve illustrate the found optimum scenario. The Central area (hypothetical) corresponds to the hypothetical case where the injection would have been maintained in the Central area during the optimization period. By promoting the injection in the North area and releasing stresses in the Central area, the found optimum strategy prevents the high seismicity rates in the Central area.

Historically, even if the injected volumes have been progressively increased in the North area, a cluster of wells with high injected volumes has been continuously operating in the Central area (see Figure 12). These persistent high injected volumes in the Central area might explain the historical surge of seismicity at Oklahoma.

4 Concluding discussion

We present a three-step constrained optimization workflow that outlines injection scenarios for maximizing injected volume under a constraint imposed by seismicity rates. It represents a physics-based predictive workflow that ensures that simulated seismicity is consistent with observed seismicity. For high-volume wastewater injection around Oklahoma, the workflow can test multiple injection scenarios to find an optimum scenario which maximizes the total volume injected while avoiding the sharp increase of seismicity as observed in 2014. This increase led to the regulatory measures and ultimately to the shut-down of the injection.

In the first step of the workflow, a forward modelling strategy is designed that honors as much as possible all available knowledge from Oklahoma. It includes information on 1) geology and flow which are used to set up flow simulations, and 2) in situ stress conditions, fault orientations, observed seismicity and prime physical processes controlling it, which are used to deploy the geomechanical and seismological analysis. Flow simulations are performed by using the open source OPM-FLOW simulator which uses historical monthly injection rates of each well at Oklahoma as input. With the approach, robust spatio-temporal pressure distributions have been computed, including complex flow interaction between the 220 wells at Oklahoma. With the simulations, the pore pressure at the nucleation depth of seismic events in the basement was accurately captured. Based on the depth distribution of the observed events, an average depth of 4 km was selected for our analysis but note here that selecting an average nucleation depth slightly deeper, e.g., 5 km, would not have affected the conclusions of our analysis. The changes in pore pressure are the main drivers of the two physical processes controlling the nucleation of induced earthquakes in Oklahoma: the direct decrease in the effective normal stress at fault due to the pore pressure increase (Langenbruch and Zoback, 2016; Norbeck and Horne, 2016; Dempsey and Riffault, 2019), and the poro-elastic loading caused by the changes in the rock volume when its pore pressure is increased (Goebel et al., 2017a; Goebel et al., 2017b; Zhai et al., 2019). The stress contribution from the poro-elastic loading has been derived from the mechanical simulator MACRIS (Candela et al., 2019; van Wees et al., 2019). MACRIS is based on a newly developed mesh-free approach, which can 1) directly use 3D pressure fields from OPM-FLOW as input for geomechanical modelling, and 2) provide high stress resolution at the interface of interest (the horizontal-plane at 4 km depth in the Oklahoma case). Coulomb stress changes have been combined with Dieterich’s theory (Dieterich, 1994) to model the spatio-temporal evolution of the seismicity rate.

In the second step, the seismicity data are assimilated in order to update the model parameters of the forward model from the first step. Data assimilation is used to assess the predictive power of the forward model by comparing simulated and historical seismicity rates. This assimilation step has often not been considered in previous studies of induced seismicity at Oklahoma injection sites. Generally, a sensitivity analysis for temporal predictions of induced seismicity is performed (e.g., Zhai et al., 2019). In the current workflow, the seismicity data has been assimilated in both space and time (see also Candela et al., 2019). It has been demonstrated that the distribution of seismicity in both space and time can be used to constrain the posterior distributions of model parameters. We showed that 1) the field-wide modelled yearly event rate and 2) the modelled spatio-temporal event densities are both successfully capturing the spatio-temporal distributions of observed events. These results confirmed the predictive power of our modelling strategy that aims to honor the ensemble of available information for the Oklahoma injection sites. Although not performed in the present study, the results suggest that the approach is suitable to forecast the potential return of the Oklahoma induced seismicity to the background rate following the arrest of all injection activities (Langenbruch et al., 2018; Zhai et al., 2019).

In the third step, an optimization approach is outlined that aims to find an optimum spatio-temporal injection scheme in order to maximize the total volume injected while keeping the seismicity rate below a cap. The cap is chosen so that the stop of the injection activities in April 2017 may be prevented as seismicity remains below a threshold value.

The present study combines both cutting-edge physics-based predictive models with a cutting-edge optimization algorithm. Despite that the modelling framework is applied to the Oklahoma case study, the primary objective was to demonstrate that complex optimization problems with two conflicting objectives involving full physics-based models for flow, geomechanics and induced seismicity can be solved. The modelling framework is based on existing workflows that were originally deployed for the optimization of well planning for hydrocarbon recovery (Chen et al., 2009; Fonseca et al., 2014, 2017), and adjusted for our specific constrained optimization problem. As such this modelling framework can thus be seen as generic and can be applied to other instances of anthropogenic subsurface activities as for example but not limited to that carbon storage and sequestration.

We showed that it likely would have been possible to avoid the dramatic rise of the rate of seismicity starting in 2014 while still injecting a total volume of fluid identical to the historical injected volume. The optimum strategy involved more uniform spatio-temporal distributions of the injection rates. More specifically the optimizer suggested to prevent the injection of additional large volumes in the Central area because already at steady-state seismicity rate due to large fluid-volumes injected before 2014.

More constraints should be added to the present approach in order to include additional key factors which have influenced the spatio-temporal historical distribution of the injection rates. As an example, the use of wells for injection should be constrained by additional parameters such as the supply of hydraulic fracturing and formation fluids from nearby hydrocarbon industry activities and notably shale gas sites as the fluid brine injected in the Arbuckle aquifer is a waste-product of production of shale gas. Accordingly, a spatio-temporal correlation between injection volumes and hydrocarbon industry operations in the area is likely. In the current optimization example, it is assumed that all wells were available at any moment in time. One way to honor the correlation between hydrocarbon industry operations and injection is to add a cost constraint in the optimization, i.e., injections via wells that were not used in historical injection should be penalized with a higher cost. This additional implementation can be achieved with the present optimization algorithm.

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions

TC: Conception and design of the work; Acquisition, analysis, and interpretation of data for the work; Drafting the work. CG: Conception and design of the work; Acquisition, analysis, and interpretation of data for the work; Drafting the work. OL: Conception and design of the work; Interpretation of data for the work; Drafting the work. JT: Drafting the work.

Acknowledgments

This work is part of the SECURE project. The data generated by the integrated approach are available on request from the corresponding author. We thank the Oklahoma Corporate Commission (OCC) and Oklahoma Geological Survey (OGS) for making publicly available respectively the disposal well data and the earthquake catalog (https://www.ou.edu/ogs).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Ader, T., Lapusta, N., Avouac, J.-P., and Ampuero, J. P. (2014). Response of rate-and-state seismogenic faults to harmonic shear-stress perturbations. Geophys. J. Int. 198 (1), 385–413. doi:10.1093/gji/ggu144

CrossRef Full Text | Google Scholar

Alt, R. C., and Zoback, M. D. (2016). In situ stress and active faulting in Oklahoma. Bull. Seismol. Soc. Am. 107 (1), 216–228. doi:10.1785/0120160156

CrossRef Full Text | Google Scholar

Baiesi, M., and Paczuski, M. (2004). Scale-free networks of earthquakes and aftershocks. Phys. Rev. E 69, 066106. doi:10.1103/physreve.69.066106

PubMed Abstract | CrossRef Full Text | Google Scholar

Barnes, J., and Hut, P. (1986). A hierarchical O(N log N) force-calculation algorithm. Nature 324, 446–449. doi:10.1038/324446a0

CrossRef Full Text | Google Scholar

Candela, T., Osinga, S., Ampuero, J. P., Wassing, B., Pluymaekers, M., Fokker, P. A., et al. (2019). Depletion-induced seismicity at the Groningen gas field: Coulomb rate-and-state models including differential compaction effect. J. Geophys. Res. Solid Earth 124, 7081–7104. doi:10.1029/2018jb016670

CrossRef Full Text | Google Scholar

Chen, Y., Oliver, D. S., and Zhang, D. (2009). Efficient ensemble-based closed-loop production optimization. SPE J. 14 (4), 634–645. doi:10.2118/112873-pa

CrossRef Full Text | Google Scholar

Dempsey, D., and Suckale, J. (2016). Collective properties of injection-induced earthquake sequences: 1. Model description and directivity bias. J. Geophys. Res. 121 (5), 3609–3637. doi:10.1002/2015JB012550

CrossRef Full Text | Google Scholar

Dempsey, D., and Riffault, J. (2019). Response of induced seismicity to injection rate reduction: Models of delay, decay, quiescence, recovery, and Oklahoma. Water Resour. Res. 55, 656–681. doi:10.1029/2018WR023587

CrossRef Full Text | Google Scholar

Dieterich, J. H., and Kilgore, B. (1996). Implications of fault constitutive properties for earthquake prediction. Proc. Natl. Acad. Sci. U. S. A. 93 (9), 3787–3794. doi:10.1073/pnas.93.9.3787

PubMed Abstract | CrossRef Full Text | Google Scholar

Dieterich, J. H. (1994). A constitutive law for rate of earthquake production and its application to earthquake clustering. J. Geophys. Res. 99, 2601–2618. doi:10.1029/93jb02581

CrossRef Full Text | Google Scholar

Ellsworth, W. L. (2013). Injection-induced earthquakes. Science 341, 1225942. doi:10.1126/science.1225942

PubMed Abstract | CrossRef Full Text | Google Scholar

Faith, J. R., Blome, C. D., Pantea, M. P., Puckette, J. O., Halihan, T., Osborn, N., et al. (2010). Three-dimensional geologic model of the Arbuckle-Simpson Aquifer, south-central Oklahoma. US Geological Survey. USGS report. Available at: https://pubs.usgs.gov/of/2010/1123/downloads/Report/OF10-1123.pdf.

Google Scholar

Firkins, M., Kolawole, F., Marfurt, K. J., and Carpenter, B. M. (2020). Attribute assisted characterization of basement faulting and the associated sedimentary sequence deformation in north-central Oklahoma. Interpretation 8, SP175–SP189. doi:10.1190/INT-2020-0053.1

CrossRef Full Text | Google Scholar

Fonseca, R. M., Stordal, A. S., Leeuwenburgh, O., Van den Hof, P. M. J., and Jansen, J. D. (2014). “Robust ensemble-based multi-objective optimization,” in ECMOR XIV-14th European conference on the mathematics of oil recovery. Sicily, Italy, September 8–11, 2014.

CrossRef Full Text | Google Scholar

Fonseca, R. R. M., Chen, B., Jansen, J. D., and Reynolds, A. (2017). A stochastic simplex approximate gradient (StoSAG) for optimization under uncertainty. Int. J. Numer. Meth. Engng. 109 (13), 1756–1776. doi:10.1002/nme.5342

CrossRef Full Text | Google Scholar

Geertsma, J. (1973). A basic theory of subsidence due to reservoir compaction: the homogeneous case. Verh. Kon. Ned. Geol. Mijnbouwkd. Genoot. 28, 43–62.

Google Scholar

Goebel, T. H. W., Weingarten, M., Chen, X., Haffner, J., and Brodsky, E. E. (2017a). The 2016 Mw5.1 Fairview, Oklahoma earthquakes: Evidence for long-range poro-elastic triggering at >40 km from fluid disposal wells. Earth Planet. Sci. Lett. 472, 50–61. doi:10.1016/j.epsl.2017.05.011

CrossRef Full Text | Google Scholar

Goebel, T. H. W., Walter, J., Murray, K., and Brodsky, E. E. (2017b). Comment on: How will induced seismicity in Oklahoma respond to decreased saltwater injection rate. Sci. Adv. 3, e1700441. doi:10.1126/sciadv.1700441

PubMed Abstract | CrossRef Full Text | Google Scholar

Hicks, A. (2011). Clustering in multidimensional spaces with applications to statistical analysis of earthquake clustering. MSc Thesis. Reno: Department of Mathematics and Statistics, University of Nevada.

Google Scholar

Hincks, T., Aspinall, W., Cooke, R., and Gernon, T. (2018). Oklahoma’s induced seismicity strongly linked to wastewater injection depth. Science 359, 1251–1255. doi:10.1126/science.aap7911

PubMed Abstract | CrossRef Full Text | Google Scholar

Holland, A. A. (2013). Optimal Fault orientations within Oklahoma. Seismol. Res. Lett. 84 (5), 876–890. doi:10.1785/0220120153

CrossRef Full Text | Google Scholar

Johann, L., Shapiro, S. A., and Dinske, C. (2018). The surge of earthquakes in Central Oklahoma has features of reservoir-induced seismicity. Sci. Rep. 8, 11505. doi:10.1038/s41598-018-29883-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, K. S. (2008). “Geologic history of Oklahoma,” in Earth sciences and mineral resources of Oklahoma (Oklahoma Geological Survey Educational Publication), 9, 3–5. 346.

Google Scholar

Keranen, K. M., Savage, H. M., Abers, G. A., and Cochran, E. S. (2013). Potentially induced earthquakes in Oklahoma, USA: Links between wastewater injection and the 2011 Mw 5.7 earthquake sequence. Geology 41 (6), 699–702. doi:10.1130/g34045.1

CrossRef Full Text | Google Scholar

Keranen, K. M., Weingarten, M., Abers, G. A., Bekins, B. A., and Ge, S. (2014). Sharp increase in central Oklahoma seismicity since 2008 induced by massive wastewater injection. Science 345, 448–451. doi:10.1126/science.1255802

PubMed Abstract | CrossRef Full Text | Google Scholar

Kolawole, F., Johnson, C., Chang, J. C., Marfurt, K. J., Lockner, D. A., Reches, Z., et al. (2019). The susceptibility of Oklahoma’s basement to seismic reactivation. Nat. Geosci. 12, 839–844. doi:10.1038/s41561-019-0440-5

CrossRef Full Text | Google Scholar

Langenbruch, C., and Zoback, M. D. (2016). How will induced seismicity in Oklahoma respond to decreased saltwater injection rates? Sci. Adv. 2, e1601542. doi:10.1126/sciadv.1601542

PubMed Abstract | CrossRef Full Text | Google Scholar

Langenbruch, C., Weingarten, M., and Zoback, M. D. (2018). Physics-based forecasting of manmade earthquake hazards in Oklahoma and Kansas. Nat. Commun. 9, 3946. doi:10.1038/s41467-018-06167-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Marone, C. (1998). Laboratory-derived friction laws and their application to seismic faulting. Annu. Rev. Earth Planet. Sci. 26, 643–696. doi:10.1146/annurev.earth.26.1.643

CrossRef Full Text | Google Scholar

McGarr, A., and Barbour, A. J. (2017). Wastewater disposal and the earthquake sequences during 2016 near Fairview, Pawnee, and Cushing, Oklahoma. Geophys. Res. Lett. 44, 9330–9336. doi:10.1002/2017gl075258

CrossRef Full Text | Google Scholar

Mindlin, R. D. (1936). Force at a point in the interior of a semi infinite Solid. Physics 7, 195–202. doi:10.1063/1.1745385

CrossRef Full Text | Google Scholar

Noda, H., Dunham, E. M., and Rice, J. R. (2009). Earthquake ruptures with thermal weakening and the operation of major faults at low overall stress levels. J. Geophys. Res. 114, B07302. doi:10.1029/2008JB006143

CrossRef Full Text | Google Scholar

Norbeck, J. H., and Horne, R. N. (2016). Evidence for a transient hydromechanical and frictional faulting response during the 2011M_w 5.6 Prague, Oklahoma earthquake sequence. J. Geophys. Res. Solid Earth 121, 8688–8705. doi:10.1002/2016JB013148

CrossRef Full Text | Google Scholar

Ogata, Y. (1998). Space-time point-process models for earthquake occurrences. Ann. Inst. Stat. Math. 50 (2), 379–402. doi:10.1023/a:1003403601725

CrossRef Full Text | Google Scholar

Okada, Y. (1992). Internal deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc. Am. 82 (2), 1018–1040. doi:10.1785/bssa0820021018

CrossRef Full Text | Google Scholar

Pei, S., Peng, Z., and Chen, X. (2018). Locations of injection-induced earthquakes in Oklahoma controlled by crustal structures. J. Geophys. Res. Solid Earth 123, 2332–2344. doi:10.1002/2017jb014983

CrossRef Full Text | Google Scholar

Qin, Y., Chen, X., Walter, J. I., Haffener, J., Trugman, D. T., Carpenter, B. M., et al. (2019). Deciphering the stress state of seismogenic faults in Oklahoma and southern Kansas based on an improved stress map. J. Geophys. Res. Solid Earth 124, 12920–12934. doi:10.1029/2019JB018377

CrossRef Full Text | Google Scholar

Rasmussen, A. F., Sandve, T. H., Bao, K., Lauser, A., Hove, J., Skaflestad, B., et al. (2021). The open porous media flow reservoir simulator. Comput. Math. Appl. 81, 159–185. doi:10.1016/j.camwa.2020.05.014

CrossRef Full Text | Google Scholar

Schmitt, S. V., Segall, P., and Dunham, E. M. (2015). Nucleation and dynamic rupture on weakly stressed faults sustained bythermal pressurization. J. Geophys. Res. Solid Earth 120 (11), 7606–7640. doi:10.1002/2015JB012322

CrossRef Full Text | Google Scholar

Segall, P., and Lu, S. (2015). Injection-induced seismicity: Poro-elastic and earthquake nucleation effects. J. Geophys. Res. Solid Earth 120, 5082–5103. doi:10.1002/2015JB012060

CrossRef Full Text | Google Scholar

van Wees, J. D., Osinga, S., Van Thienen-Visser, K., and Fokker, P. A. (2018). Reservoir creep and induced seismicity: inferences from geomechanical modeling of gas depletion in the groningen field. Geophys. J. Int. 212, 1487–1497. doi:10.1093/gji/ggx452

CrossRef Full Text | Google Scholar

van Wees, J. D., Pluymaekers, M., Osinga, S., Fokker, P. A., Van Thienen-Visser, K., Orlic, B., et al. (2019). 3-D mechanical analysis of complex reservoirs: A novel mesh-free approach. Geophys. J. Int. 219, 1118–1130. doi:10.1093/gji/ggz352

CrossRef Full Text | Google Scholar

Walsh, F. R., and Zoback, M. D. (2015). Oklahoma’s recent earthquakes and saltwater disposal. Sci. Adv. 1, e1500195. doi:10.1126/sciadv.1500195

PubMed Abstract | CrossRef Full Text | Google Scholar

Weingarten, M., Ge, S., Godt, J. W., Bekins, B. A., and Rubinstein, J. L. (2015). High-rate injection is associated with the increase in U.S. mid-continent seismicity. Science 348, 1336–1340. doi:10.1126/science.aab1345

PubMed Abstract | CrossRef Full Text | Google Scholar

Zaliapin, I., and Ben-Zion, Y. (2013). Earthquake clusters in southern California, I: Identification and stability. J. Geophys. Res. Solid Earth 118, 2847–2864. doi:10.1002/jgrb.50179

CrossRef Full Text | Google Scholar

Zaliapin, I., and Ben-Zion, Y. (2016). Discriminating characteristics of tectonic and human-induced seismicity. Bull. Seismol. Soc. Am. 106 (3), 846–859. doi:10.1785/0120150211

CrossRef Full Text | Google Scholar

Zaliapin, I., Gabrielov, A., Keilis-Borok, V., and Wong, H. (2008). Clustering analysis of seismicity and aftershock identification. Phys. Rev. Lett. 101, 018501. doi:10.1103/PhysRevLett.101.018501

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhai, G., Shirzaei, M., Manga, M., and Chen, X. (2019). Pore-pressure diffusion, enhanced by poro-elastic stresses, controls induced seismicity in Oklahoma. Proc. Natl. Acad. Sci. U. S. A. 116 (33), 16228–16233. doi:10.1073/pnas.1819225116

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhuang, J., and Touati, S. (2015). “Stochastic simulation of earthquake catalogs,” in Community online resource for statistical seismicity analysis. Available at: http://www.corssa.org. doi:10.5078/corssa-43806322

CrossRef Full Text | Google Scholar

Zhuang, J., Harte, D., Werner, M. J., Hainzl, S., and Zhou, S. (2012). “Basic models of seismicity: temporal models,” in Community online resource for statistical seismicity analysis. doi:10.5078/corssa-79905851

CrossRef Full Text | Google Scholar