# Response of induced seismicity to injection rate reduction: models of delay, decay, quiescence, recovery and Oklahoma

This Jupyter Notebook provides an executable and interactive interface to selected models from Dempsey and Riffault *"Response of induced seismicity to injection rate reduction: models of delay, decay, quiescence, recovery and Oklahoma"*. We provide an abridged explanation of the model, sufficient to introduce the equations and figures. Further discussion and justification of the model assumptions can be found in the full version of the article.

The models are implemented in Python. The simplest way to get them running is to install the [Anaconda Python](https://www.anaconda.com/download/) distribution, Version 2.X. To run models using the FEHM reservoir simulator, ensure you have downloaded a copy of [FEHM](https://github.com/lanl/FEHM) and installed [PyFEHM](http://pyfehm.lanl.gov/).

## Pressure around a point source injector

Triggering of induced seismicity is controlled by the build-up of fluid pressure underground. In many instances, this is achieved by injection into geologic formations whose horizontal extent greatly exceeds their vertical thickness, and that are bounded above and below by much less permeable formations. For this geometry, flow is primarily radial (Fig. 1B). Solving for mass balance under conditions of Darcy flow and homogeneous material properties gives the pressure change at radius, $r'$, and time, $t'$, due to injection at constant mass rate, $q'$, as

\begin{equation}
p'-p'_0 = \frac{1}{II'}q' W\left( \frac{r'^2}{4D' t'} \right),\quad II'=\frac{4\pi k'h'\rho'}{\eta'},\quad D'=\frac{k'}{\eta'(\phi \beta'_{fl}+\beta_{aq}')},\quad W(x) = -\int\limits_{-x}^{\infty} \frac{e^{-\xi}}{\xi} d\xi,
\end{equation}

where $p'_0$ is the initial pressure, $II'$ is the injectivity index of the well, $\eta'$, $\rho'$ and $\beta'_{fl}$ are, respectively, fluid (dynamic) viscosity, density and compressibility, $D'$, $k'$, $h'$, and $\beta'_{aq}$ are, respectively, aquifer diffusivity, permeability, thickness, and compressibility, and $W(x)$ is the exponential integral. For small values of $x$, the exponential integral can be approximated $W(x)\approx 0.5772+\ln(x)$. A key feature of the equations above is that, for all $r'$, pressure increases monotonically and, eventually, logarithmically. This simple model must eventually break down (it predicts infinite pressure), typically when $p'$ becomes so large as to overwhelm - possibly by brittle failure - the containment effect of the upper and lower bounding formations.

We consider the particular situation that, after a period of injection, $t'_{cut}$, the injection rate is instantaneously reduced by a fraction, $f$, to a new constant rate, $(1-f)q'$. Pressure evolution is then obtained by superposition of the solutions above for two injectors, one beginning at $t'=0$ at rate $q'$, and the second beginning at $t'=t'_{cut}$ at the negative rate $-fq'$ (i.e., a producing well):

\begin{equation}
p'-p'_0 = \frac{1}{II}q' \left[W\left( \frac{r'^2}{4D' t'} \right) - f\,W\left( \frac{r'^2}{4D' (t'-t'_{cut})} \right)\right]\quad\quad(1).
\end{equation}

It is convenient to rewrite the pressure equation in a dimensionless form using the dimensionless quantities $t = t'/t'_{cut}$, $r = r'/2\sqrt{D' t'_{cut}}$ and $p = (p'-p'_0)II'/q'$. Here, the diffusion length scale, $r'_D$, has been used to dimensionalize distance, with $t'_{cut}$ the characteristic time for the system. The resulting pressure equation is

\begin{equation}
p=W\left( \frac{r^2}{t} \right)-fW\left( \frac{r^2}{t-1} \right),
\end{equation}

where $f=0$ for $t<1$. Throughout this article, we shall adopt the convention that primed variables are dimensional, and unprimed variables are dimensionless.

***Execute the cell below (select it $\rightarrow$ Ctrl+Enter) to plot pressure as a function of distance and time. Use slider bars and dropdown to change the plot.***


```python
%matplotlib inline
from plotting import *

# pressure plotting function
def pressure(show='distance', advance=1, f=0.4):  
    iso = InducedSeismicity(f, dim=2)              # model object
    fig, ax = plt.subplots(1)                      # create figure
    draw_pressure(fig,ax,iso,show,advance)         # add plot
    plt.show()                                     # display

# interactive plotting
interact(pressure,  show = ['distance','time','pressure'], advance = (1,9,1), f=(0.1,0.9,0.1));
```

## Pressure around a distributed circular injector

A useful generalization of the Theis equation is obtained for the situation where mass injection is evenly distributed inside the radial area $r'<r'_{inj}$. Introduction of a second length scale provides a choice of dimensionalization and, for these models (including those incorporating a seismogenic basement), we have opted to define the dimensionless radius, $r=r'/r'_{inj}$. 

The solution to this problem is given by *Carslaw and Jaeger* [1959] $\S$13.8, Eqs. (13)-(14)
\begin{equation}
p_1= \frac{16}{\pi^2}\int\limits_0^\infty \frac{(1-f-e^{-(r_Du)^2t}+fe^{-(r_Du)^2(t-1)})J_0(ru)J_1(u)}{u^4\left[J_1(u)Y_0(u) - J_0(u)Y_1(u)\right]^2}\,du
\end{equation}

\begin{equation}
p_2 = \frac{8}{\pi}\int\limits_0^\infty \frac{(1-f-e^{-(r_Du)^2t}+fe^{-(r_Du)^2(t-1)})J_0(ru)J_1(u)}{u^3\left[J_1(u)Y_0(u) - J_0(u)Y_1(u)\right]}\,du 
\end{equation}

where $p_1$ and $p_2$ are, respectively, pressure inside and outside the injection radius, $r_D$ is the dimensionless diffusion length, equivalent to the ratio of diffusion and injection length scales (i.e., $r_D=\sqrt{D' t'_{cut}}/r'_{inj}$), $J_i$ and $Y_i$ are, respectively, $i^\text{th}$ order Bessel functions of the first and second kind, and, as in the previous section, $f=0$ for $t<1$. 

***Execute the cell below to plot pressure as a function of distance. Use slider bars to change the plot.***


```python
# pressure plotting function
def pressure_distributed(loga=-1, f=0.4, advance=1):  
    fig, ax = plt.subplots(1)                           # create figure
    draw_pressure_distributed(fig,ax,loga,f,advance)    # add plot
    plt.show()                                          # display

# interactive plotting
interact(pressure_distributed,  loga=(-1,1,0.5), f=(0.1,0.9,0.1), advance = (1,9,1));
```

## Steady-state stressing model of induced seismicity


At steady-state, the seismicity rate, $n_{ss}$, implied by a 0-dimensional rate-state model [[*Dieterich*, 1994](https://agupubs.onlinelibrary.wiley.com/doi/10.1029/93JB02581);[*Segall and Lu*, 2015](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2015JB012060)] is proportional to the pressure loading, $dp/dt$,
\begin{equation}
n_{ss} \propto \frac{dp}{dt},
\end{equation}

In this study, we parameterise fault proximity to failure using a Mohr-Coulomb criterion coupled to a stressing rate model. The principle steps we follow to compute a global seismicity rate, $N(t)$, for a system with arbitrary pressure evolution, $p(r,t)$, are:

1. Assume an infinite, homogeneous medium through which identically oriented and pre-stressed faults are uniformly distributed, as in [*Segall and Lu* [2015]](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2015JB012060). 
2. Define a critical pressure rise, $p_c$. The local seismicity rate, $n(r,t)$, is zero until until $p_c$ is exceeded.
3. Once exceeded, $n(r,t)$ is proportional to the local pressure loading, $\partial p/\partial t$, as given by the steady-state stressing model. Proportionality is assumed constant in time. 
4. Compute the global seismicity rate, $N(t)$, as the integral over space of the local seismicity rate. For a point injection source, this is evaluated semi-analytically. For a distributed injection configuration, $p$ is determined by a reservoir simulator, $\partial p/\partial t$ is evaluated by finite differences, and the integral over space is computed numerically using the Trapezium rule.
5. If pressure should locally decrease in response to an injection rate change, then the local seismicity rate is immediately zero. Pressure must exceed the previous maximum value at that location before seismicity can restart, at which time it is once again proportional to the local pressure loading. The retreat from and then return to failure is sometimes referred to as the Kaiser effect [[*Baisch and Harjes*,2003](https://onlinelibrary.wiley.com/doi/abs/10.1046/j.1365-246X.2003.01837.x)].

## Induced seismicity around a point source injector

Following points 3 and 4 above and using the pressure results developed earlier for point injection into a 2D geometry, the global seismicity rate prior to an injection reduction is expressed

\begin{equation}
N = \kappa \int\limits_0^{r_c}2\pi r\, \frac{\partial p}{\partial t} dr = \pi\kappa\left(1-e^{-r_c^2/t}\right),
\end{equation}

where $r_c(t)$ is the distance from the well at which fluid pressure has risen to the critical pressure, i.e., $r_c = \sqrt{t W^{-1}(p_c)}$ ,and $r_c\rightarrow\infty$ as $p_c\rightarrow 0$. Therefore, at constant injection rate, the seismicity rate is constant for any $p_c$. 

We explicitly introduce a constant of proportionality, $\kappa$, which modulates the amount of seismicity per unit pressure increase over the integrated area. Although $\kappa$ is not defined in terms of other properties of the crust or injected fluid, it nevertheless aggregates the modulating effect of fracture density, seismogenic fault depth, degree of aseismic creep, etc., on the seismic productivity of fractured rock under fluid load.

Once injection is reduced, pressure starts to decline at some locations and seismicity temporarily ceases in these regions. For $t>1$, an adjustment is made to the limits of the seismicity integral
\begin{equation}
N = \kappa \int\limits_{r_{in}}^{r_c} 2\pi r\,\frac{\partial p}{\partial t} dr = -\pi\kappa\left[ e^{-r^2/t} -fe^{-r^2/(t-1)}\right]_{r_{in}}^{r_{c}},
\end{equation}

to exclude those regions in which $\partial p/\partial t<0$ or where fluid pressure has yet to exceed a previous maximum (Kaiser effect). Finally, the seismicity rate is normalized by its value immediately prior to the rate reduction (i.e., $N=1$ at $t=1$), which drops the term $\kappa$ from subsequent expressions.

Two limit schemes are identified: 

1. For $p_c>p^*$ then, although seismicity may temporarily cease at some locations, there are others that generate seismicity at all times, i.e., the global seismicity rate is always non-zero. In this case, $r_{in}=g(t)$ for $t<t^*$, and $r_{in}=h(t)$ for $t>t^*$.
2. For $p_c<p^*$, there will be a finite period of time for which seismicity has ceased everywhere around the well, only to start up again later. We delimit the quiescent period $[t_g, t_h]$, with the boundaries respectively determined by numerical inversion of $p_c=p(g(t_g),t_g)$ and $p_c=p(h(t_h),t_h)$. Then, as in the first case, $r_{in}=g(t)$ for $t<t_g$, and $r_{in}=h(t)$ for $t>t_h$, and $N=0$ for $t_g<t<t_h$.  

***Execute the cell below and use the slider bar to explore how the two limit schemes arise.***


```python
# seismicity plotting function
def seismicity(pc=0.1):  
    fig, axs = plt.subplots(2,1)                      # create figure
    iso = InducedSeismicity(f=0.4, dim=2, pcrit=pc)   # model object
    draw_seismicity(fig,axs,iso)                      # add plot
    plt.show()                                        # display

# interactive plotting
interact(seismicity,  pc = (0.,2.,0.1));
```

## Sensitivity to rate cut and critical pressure

For the case of injection at a point, our model predicts that the seismicity rate should be constant prior to any changes of the injection rate. Following the rate reduction, a new (constant) seismicity rate will eventually be established following a period of equilibration. In practical situations, deviations from a stable equilibrium seismicity rate are anticipated, attributable to breakdowns in the assumption of evenly distributed and pre-stressed fractures, or homogeneous material properties. 

For the special case of faults that are right on the limit of failure ($p_c$=0), imposing a rate reduction, $f$, results in a proportionate decrease in the new equilibrium seismicity rate. This is in agreement with the result from [*Shapiro and Dinske* [2009]](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2008JB006145) that seismicity and injection rate are proportional. In the case of non-critically stressed faults ($p_c>0$), the reduction of the equilibrium seismicity rate is greater than the proportional reduction in injection.

***Execute the cell below and use the sliders to see how post rate-reduction seismicity varies with $f$ and $p_c$.*** 


```python
# seismicity plotting function
def sensitivity(f=0.4, pc=0.):  
    fig, ax = plt.subplots(1,1)                      # create figure
    iso = InducedSeismicity(f=f, dim=2, pcrit=pc)    # model object
    draw_sensitivity(fig,ax,iso)                     # add plot
    plt.show()                                       # display

# interactive plotting
interact(sensitivity,  f=(0.1,1.0,0.1), pc = (0.,2.5,0.1));
```

## Oklahoma model

To understand physical drivers behind the recent decline of seismicity in Oklahoma, we have developed a simplified model of injection, pressure build-up, and induced seismicity for the  Western Oklahoma (WO) region.  We adapted our model of 2D distributed injection overlying seismogenic basement to account for the varying injection rate since 2004. Injected mass is evenly distributed inside an circular area with radius 50$\,$km. We approximate the Arbuckle formation as a 500$\,$m thick layer of unknown isotropic permeability, $k_a'$, 10% porosity, and matrix compressibility 2$\times$10$^{-11}$Pa$^{-1}$. The basement has unknown vertical permeability, $k_b'$, and is assigned negligible horizontal permeability, consistent with a conceptual understanding of vertically oriented faults as the primary conduits conducting overpressure to seismogenic depths. The basement is assigned zero compressibility, and 1% porosity. We compute seismicity using basement pressures weighted by an envelope that follows the depth distribution of seismicity in Oklahoma.

Seismicity is computed assuming a steady-state stressing rate model and an unknown constant critical pressure, $p_c$. The scaling parameter $\kappa$ is unknown, to be estimated during calibration (note $\kappa$ values here differ from the article text by a factor of $2\pi$). Finally, a background seismicity rate is added to the simulated triggering rate, equal to 1.16$\,$M$\geq$3 event per year for the entire state of Oklahoma, scaled appropriately to the area of WO. Further details of the model can be found in the main text.

Because our model takes too long to execute for smooth use of the sliders (around eight seconds), we pre-executed 7203 forward runs covering all combinations of seven different values per parameter and three forecast scenarios. Results are stored in a .csv file and read in the cell below. 

Overlaid on profiles of the seismicity rate are 5 to 95% Poisson uncertainty bounds (shaded interval), i.e., in a given month, a random Poisson process with modeled rate parameter of 10 events per month has a 5% likelihood of generating less than 5 events 5% likelihood of generating more than 15 events.

***Execute the cell below and use the sliders to calibrate the model to the observations.*** 


```python
# load dataset from csv file
dataset = np.genfromtxt('LLK.csv', delimiter = ',')       # load saved dataset
N_models = int(dataset[0,0])                              # number of simulations pre-runned

# model plotting function
def oklahoma(ka=-11.7, kb=-15.4, pc=0.045, kappa=9.):  # default at maximum LLK
    fig, ax = plt.subplots(1,1)                           # create figure
    draw_oklahoma(fig, ax, [ka, kb, pc, kappa], dataset[:,:N_models+1]) # draw plot
    plt.show()                                            # display

# interactive plotting
interact(oklahoma,  
        ka=(-12.6, -10.8, 0.3),              # range of Arbuckle formation permeability values
        kb=(-15.55, -15.25, 0.05),              # range of basement vertical permeability values
        pc=(0.015, 0.075, 0.01),              # range of critical pressure values
        kappa=(3., 15., 2.),  # range of scaling parameter values
        );
```

22912eebc57775737c2d8637d564e2217989f681
ipynb
Jupyter Notebook
induced_seismicity_models.ipynb
ddempsey/dempsey2018_agu
c0ed1d471004353a4f525776526bf7f8b256f384
MIT
2018-12-12T13:55:21.000Z
2021-01-13T02:57:15.000Z
induced_seismicity_models.ipynb
ddempsey/dempsey2018_agu
c0ed1d471004353a4f525776526bf7f8b256f384
MIT
induced_seismicity_models.ipynb
ddempsey/dempsey2018_agu
c0ed1d471004353a4f525776526bf7f8b256f384
MIT
Qwen/Qwen-72B
1. YES
2. YES