Geo-temporal modelling: induced earthquakes in Oklahoma

For 2005-08 this catalogue lists one to three magnitude-3+ earthquakes a year in and around Oklahoma. For 2015 it lists more than nine hundred. The surge is widely attributed to the deep disposal of wastewater from oil and gas production, and it receded after regulators limited injection volumes in 2015-16. The earthquakes did not rise and fall everywhere at once: the active area moved.

This notebook asks three questions of the catalogue alone:

1. Where and when was the earthquake rate elevated?
2. How did the hot-spot move?
3. What does the model forecast for a held-out final period (2018-19), and how much of that forecast deserves to be believed?

It builds on E05: you should already know what an HSGP is, what m and c do, and why centred and non-centred coefficients behave differently. Here the GP becomes two-dimensional, gets a time axis, and sits behind a count likelihood.

#1 · Question and data

Seven events in the four years 2005-08, then 20, 41, 63 ... 914 in 2015, and back down to 60 in 2019. The largest events give a first idea of the geography:

#2 · From a catalogue to a space x time count array

A point process model needs three pieces of data engineering, and each involves a choice that ends up in the results, so we make them in the open.

Project to kilometres. A GP lengthscale in degrees means different distances east-west and north-south. With no GIS stack in this environment we use an equirectangular projection about a reference point $(\lambda_0, \varphi_0)$,

$$x = R\,\cos\varphi_0\,(\lambda - \lambda_0), \qquad y = R\,(\varphi - \varphi_0),$$

with longitudes and latitudes in radians and $R = 6371$ km. It ignores the convergence of the meridians; the cell below prints how much that matters over our window.

About 1% of distortion at the northern and southern edges: irrelevant next to 25 km cells.

Choose a window and a grid. The catalogue covers a box from Texas to Kansas, but most of it is empty. We model a 300 km x 225 km window around the active area with 25 km cells (12 x 9 = 108 cells), and quarters as the time unit. Smaller cells would give prettier maps and mostly-zero data; bigger cells would blur hot-spots that are only a few tens of km across. The clock starts in 2009, when the rate first lifts off.

Hold out the end. 2009-2017 (36 quarters) is training data; 2018-2019 (8 quarters) is never shown to any model and is used once, in section 10.

From here on everything carries the dims ("cell", "period"). An xarray DataArray makes the bookkeeping trivial: counts.sum("cell") is the time series, counts.groupby("year") the yearly maps.

Now some plotting helpers. We have no basemap, so a handful of towns provide orientation. Their coordinates agree with the catalogue's own place strings ("3 km W of Cushing").

The colours already tell the story of question 2: the early (dark) events sit in the south-east around Prague and Oklahoma City, the later (green-yellow) ones further north and west. The sharp upper edge is not geology: the catalogue query stops at latitude 37.1, just north of the Kansas border, and our window stops with it. Whatever the model says about the northernmost row of cells is about a truncated hot-spot.

The same thing as the model will see it - counts per cell, for six of the nine training years, on a logarithmic colour scale (white cells are empty):

In 2010 there is a small patch east of Oklahoma City; 2011 is dominated by the magnitude-5.7 Prague sequence in the south-east; 2013 adds a separate patch near the Kansas border; by 2014-15 activity fills a band from Oklahoma City to the border, with an arm reaching west past Cherokee; by 2017 the map is dimmer nearly everywhere. A model whose spatial pattern cannot change over time has no chance.

#3 · The model family: a log-Gaussian Cox process

A Poisson process with intensity $\lambda(s, t)$ (events per km² per unit time) says that counts in disjoint regions of space-time are independent Poisson variables whose means are the integrals of $\lambda$. A Cox process makes $\lambda$ itself random, and a log-Gaussian Cox process gives $\log\lambda$ a Gaussian-process prior. On our grid, with $\lambda$ treated as constant within a cell and a quarter,

$$y_{ct} \sim \text{Poisson}(\lambda_{ct}), \qquad \log\lambda_{ct} = \alpha + g_t + f(s_c) + h(s_c, t),$$

where $g$ is a temporal component, $f$ a smooth spatial field evaluated at the cell centre $s_c$, and $h$ a space-time interaction. The grid is a quadrature rule, not part of the model: $f$ is a function on the plane. We climb a ladder:

#The spatial field: a 2-D HSGP

Everything from E05 carries over; the basis functions become products of sines, $\phi_{jk}(x, y) = \phi_j(x)\,\phi_k(y)$, and m=[m_x, m_y] sets how many per axis. The number of basis functions is the product $m_x m_y$, which is why HSGP stops being attractive beyond two or three dimensions.

Lengthscale prior, in km. The geography brackets it. The blobs on the maps above are one to three cells across, which points at lengthscales of a few tens of km. Much below the 25 km cell width a field is indistinguishable from independent cell effects, and a lengthscale above 100 km - a third of the window - would make the whole active area one smooth hump. As in E05 we ask PreliZ for the maximum-entropy inverse-gamma with 90% of its mass between 15 and 100 km, and print what the lower end means on this grid:

At $\ell = 15$ km neighbouring cells are still correlated, but only just: the grid has little to say about anything shorter, and the inverse-gamma's thin left tail keeps the sampler out of that region.

Choosing m and c. E05's heuristic approx_hsgp_hyperparams is built for one dimension, where basis functions are cheap. Here we are under a hard constraint: in rung 3 every basis function gets one coefficient per year, so the parameter count is $m_x m_y \times 11$. So we reason directly:

• The cell centres are 25 km apart, so the grid cannot see wavelengths below 50 km. Basis function $j$ has wavelength $4L/j$; anything much finer than 50 km is wasted.
• c sets the boundary $L = c \times$ half-range. Too small and the field is pinched to zero at the window's edge - right where the northern hot-spot is.

The table compares the exact Matern-5/2 covariance between cell centres with its HSGP approximation: the lowest prior variance of any cell (it should be 1; low values mean a pinched edge) and the largest error anywhere in the covariance matrix.

• m=[10, 8], c=1.3 is cheap and wrong everywhere: the edge variance is already down to 0.70 at $\ell = 35$ km.
• m=[20, 16], c=2.0 is accurate up to $\ell = 50$ km and tolerable at 100 km, but needs 320 basis functions: 3520 coefficients in rung 3.
• m=[14, 10], c=1.5 (140 basis functions, finest wavelength 59 km - about what the grid can resolve) is accurate for $\ell \le 35$ km, acceptable at 50 km, and poor at 100 km.

We take the middle option, and with it an obligation: if the posterior of $\ell$ ends up much above 50 km, the approximation is not the model we wrote down and we must come back and pay for a larger c and m. The maps suggest it will not. This is the honest way to run an HSGP on a budget: spend basis functions where the posterior is, then check that it is where you thought.

#4 · Prior predictive: LGCP priors explode

Here is the textbook LGCP of rung 2 with priors that would pass for "weakly informative" in a regression: Normal(0, 3) on the intercept, HalfNormal(3) on the field's amplitude, HalfNormal(1) on the quarterly random-walk step. What do they claim about the number of earthquakes per quarter in the window?

We summarise the expected count rather than simulated counts for a practical reason: NumPy's Poisson sampler raises an error for the means in the upper tail. A third of the prior draws say "more than ten thousand M3+ earthquakes per quarter in central Oklahoma", the median is about a thousand, and the tail runs to numbers with no physical meaning. Three things conspire:

• Everything is added on the log scale and then exponentiated, so prior standard deviations of 3 + 3 + a random walk become factors of thousands.
• The total is a sum of exponentials of the field. For a Gaussian $f$ with standard deviation $\eta$, $E[e^f] = e^{\eta^2/2}$: at $\eta = 3$ that is a factor of 90 from the field alone, and the sum is dominated by whichever cell drew the largest value.
• A random walk's variance grows with time: after 36 quarters the log-rate has wandered $\sqrt{36} = 6$ steps' worth.

#The fix: separate how many from where

The second point is structural, and tightening $\eta$ is the wrong cure: the maps show intensities that differ by orders of magnitude between cells, so $\eta$ really is large. Instead use a classical fact about Poisson counts: independent $y_{ct} \sim \text{Poisson}(\lambda_{ct})$ are equivalent to a Poisson total $N_t \sim \text{Poisson}(\Lambda_t)$ that is allocated to cells with multinomial probabilities $\pi_{ct} = \lambda_{ct} / \Lambda_t$. So we parametrise

$$\log\lambda_{ct} = \underbrace{\alpha + g_t}_{\log \Lambda_t:\ \text{how many}} + \underbrace{\log \text{softmax}_c\big(f_t(s_c)\big)}_{\log \pi_{ct}:\ \text{where}}.$$

It is still an LGCP - the log-intensity is a Gaussian field plus an offset - but now $e^{\alpha + g_t}$ is the expected number of events in the window, whatever the field does, and the field's amplitude only controls how concentrated the map is. Priors become easy to state:

• $\alpha \sim \text{Normal}(\log \bar N, 1.5)$ with $\bar N$ the mean quarterly total: the typical quarter is within a factor of $e^{\pm 3}$ of it. (Using the data to centre a weak prior is a convenience; the width is what matters.)
• $g$: a random walk with sigma_g ~ HalfNormal(0.5) per quarter, centred to average zero over the training period so that $\alpha$ keeps its meaning. The catalogue went from a handful to hundreds of events per quarter, so large drifts must stay possible.
• Field amplitude eta0 ~ HalfNormal(2): relative intensities between cells of up to $e^{\pm 4}$. Hot cells have hundreds of events and most have none.

This parametrisation will pay off a second time when we forecast (section 10).

#One builder for the whole ladder

The four rungs share almost everything, so one function builds them all. Details:

• Centred or non-centred? Both, for the reasons given in E02 and E05. Each quarter's total is observed directly, so the steps of $g$ are strongly informed: centred, step ~ Normal(0, sigma_g). Most of the 140 x 11 field coefficients are weakly informed (most cells are empty most of the time): non-centred, z ~ Normal(0, 1) scaled afterwards.
• The future is inside the model. period and era run to the end of 2019, but the likelihood only sees the training slice. The innovations of 2018-19 have no data attached, so their posterior is their prior - and every posterior draw carries its own simulated future, consistent with that draw's hyperparameters. Forecasting will need no extra machinery.
• gp.prior_linearized (mentioned in E05) returns the basis phi and the spectral weights sqrt_psd, so that we can attach our own coefficients. The cell centres never change, so they go in as a constant rather than pm.Data, and there is nothing to freeze.
• In rung 3 the field of year $e$ is $f_e = f_{e-1} + \delta_e$, where each increment $\delta_e$ is itself a draw from the spatial GP with amplitude eta_d: a random walk in function space. In the basis that is simply a random walk on each coefficient.
• A random walk needs an anchor: the one year whose map gets the prior $\text{GP}(0, \eta_0^2 k_\ell)$ instead of being a sum of increments. The first year is the habitual choice and a poor one here. 2009 has 20 events, so its map is barely identified, and every well-observed year would be a sum of five to seven increments, all of which the data then correlate. We anchor at 2015, the best-observed year, and let increments run backwards to 2009 and forwards to 2019. It is the same random walk - independent, symmetric increments - with the weakly informative marginal placed where the data are; eta0 becomes "the amplitude of the 2015 map". (Set ANCHOR_YEAR = 2009 to see what it buys: the posterior is essentially unchanged and sampling is markedly slower.)

Prior predictive on the count scale, for the most flexible rung. Now that totals are bounded we can simulate actual counts:

The new prior is still broad - it has to be, for a rate that really did rise more than a hundredfold - but it lives on the right planet: quarterly totals from a couple to about a thousand, maps in which most cell-quarters are empty, and a largest count somewhere between a handful and a few thousand. The upper tail remains generous, but the observed values sit inside the central 90% of every row, and nothing is off by ten orders of magnitude.

#5 · Rung 1: time only

Every cell gets the same share, $1/108$. The helper below samples a rung and reports the diagnostics we will want each time.

No divergences, r_hat of 1.00. To look at predictions we need the expected counts $\mu_{ct}$ for all 44 quarters and simulated counts around them. Because the future innovations are in the posterior, one helper serves both the in-sample replicates and the forecast. It is written with xarray so that the era-to-quarter lookup is done by label, and it thins to 1000 draws to keep the arrays small.

The temporal check passes almost by construction: a random walk with one free value per quarter can follow any series of totals (the axis is linear below 10 and logarithmic above). Note the spike in late 2011 - the Prague sequence, a mainshock with its aftershocks. The model can only read that as a sudden rise and fall of the background rate.

The spatial check fails completely. Spreading events evenly over 108 cells predicts far fewer empty cells than the 86% observed. On to space.

#6 · Rung 2: separable space + time

Fast and clean. eta0 is large - the log-intensity varies across the map with a standard deviation approaching 4, i.e. the busiest places are thousands of times more active than the quietest - and the lengthscale is short, about 30 km: not much more than one cell. Both are as the maps suggested, and $\ell$ is in the range where our m and c are accurate.

Is one fixed map enough? Compare observed yearly counts with the model's expectation:

The separable model paints the same picture every year and only turns the brightness up and down. It puts the 2011 Prague sequence on a faint copy of the 2015 map, predicts activity in the north-west years before there was any, and keeps the south-east glowing in 2017. The time-averaged map is a fair summary of nine years and a poor description of any one of them. That is what "no interaction" means.

#7 · Rung 3: letting the hot-spot move

We need $f$ to depend on time. Options, and why we choose the one we do:

Yearly rather than quarterly eras keep the parameter count down by a factor of four; the quarter-to-quarter changes in the total are still carried by $g_t$. The price is that the map can only change on 1 January.

We fit it first with the Poisson likelihood of the textbook LGCP.

1540 field coefficients, a few seconds of sampling, no divergences and r_hat at 1.00-1.01. The GP hyperparameters ell and eta0 mix most slowly (bulk ESS of a few hundred); that is enough for a stepping stone, and we will run the final model for longer.

eta_d says the log-intensity at a typical location changes by about 1.6 per year - a factor of five up or down - and that is in the upper tail of its HalfNormal(1) prior: the data insist. The map is not drifting, it is being redrawn.

#Is Poisson enough?

Conditional on $\lambda$, an LGCP assumes events are independent. Earthquakes are not: they trigger aftershocks, so events arrive in bursts that are tighter in space and time than any smooth field. Check two statistics the Poisson likelihood constrains only indirectly - the share of empty cell-quarters and the largest single count:

Each title gives $p$, the share of replicated datasets in which the statistic is at least as large as observed. Both checks fail, in the same direction. Replicated datasets have too few empty cells and hardly ever a count as large as the observed 86: the real events are more clumped than a Poisson process with this smooth intensity. The field has done what it could to absorb the clumps - remember how rough and volatile it came out - but a smooth surface cannot produce a burst in one cell in one quarter.

#Rung 4: negative binomial counts

$y_{ct} \sim \text{NegBinomial}(\mu_{ct}, \phi)$ has variance $\mu + \mu^2/\phi$. We put the prior on $1/\phi$ (inv_phi ~ Exponential(1)), so that the Poisson model is the point $1/\phi = 0$ rather than $\phi = \infty$. Think of it as a crude stand-in for aftershock clustering: it says bursts happen, without saying when or where.

No divergences, r_hat of 1.00-1.01 for every hyperparameter and for the whole field. This is the stickiest fit of the notebook: ell and eta0 trade off against each other and against all 1540 coefficients at once (a longer lengthscale with a larger amplitude gives a similar field), so their bulk ESS is 700-800 out of 12,000 draws, against thousands for the rest. That is why this rung gets 3000 draws per chain.

The rank plot tells the same story. (Its envelope assumes independent draws, so we thin the chains first; unthinned, every slowly mixing parameter is flagged however long you run.) Four panels stay inside the 99% envelope. ell does not quite: one chain spent less time at short lengthscales than the others and touches the edge (black marks). With r_hat at 1.00 that is slow mixing rather than chains disagreeing about where the posterior is - adequate for the maps and intervals below, but run longer if a hyperparameter is itself what you want to report.

inv_phi is far from zero: the data want the extra variance, and both predictive checks now pass. Compare the two columns of the last table for what that did to the rest of the model. With bursts handled by the likelihood, the field relaxes: a longer lengthscale (about 36 km instead of 27) and smaller yearly increments, and the random walk for the total takes smaller steps (sigma_g) because a quarter like 2011 Q4 no longer has to be explained by the background rate alone. An overdispersed likelihood is a regulariser for the latent field. The posterior of $\ell$ stays well below 50 km, where the m=[14, 10], c=1.5 approximation was shown to be accurate, so the obligation from section 3 is discharged.

#8 · Answers: where, when, and how it moved

The posterior mean of the expected number of events per cell per year, with the actual epicentres on top (blue dots):

Where and when: a faint patch east of Oklahoma City in 2009-10; the Prague sequence in 2011-12; in 2013 a new centre north-east of Oklahoma City and a first, separate patch near the Kansas border; in 2014-15 both grow and merge into a band from Oklahoma City to the border, with an arm reaching west past Cherokee; from 2016 the whole map dims. The brightest cell of 2016-17 lies north-west of Fairview, which is where the catalogue places the magnitude-5.1 event of February 2016. White cells have fewer than 0.3 expected events a year.

A mean is half an answer. The map of uncertainty - here the posterior standard deviation of the log expected count - is just as informative:

The pattern is the mirror image of the intensity map. Where there are events the field is pinned down to a fraction of a log unit; in empty regions all the data say is "low", and how low is left to the prior - an sd of 2 on the log scale means "somewhere between very few and almost none". That is the right kind of ignorance, and it is worth showing to anyone tempted to read the blank cells of the mean map as a precise zero. In 2010 - a few dozen events, and five random-walk steps away from the anchor year - everything outside the patch near Oklahoma City is uncertain, much of it beyond the 2.6 at which the colour scale is capped.

How did the hot-spot move? Summarise each year's map by its centre of mass, $\sum_c \pi_{c,e}\, s_c$, a derived quantity with a full posterior. The observed counterpart is the centroid of that year's events (by cell centre, to compare like with like).

The interaction model follows the observed centroid year by year, with wider bands in the sparse early years and tight ones at the peak. Between the Prague year and the 2015 peak the centre of activity moved well over 100 km to the north-west - some 65 km west and 100 km north, as printed above - and then barely moved while the rate fell. The separable model, by construction, has a single centroid for all years (its band is flat). To the right of the dotted line the bands are forecasts; they open up quickly, which is the subject of section 10.

#Spatial posterior predictive check

The same observed-versus-expected maps as in section 6, now with rung 4 added:

Rung 4 reproduces each year's map, including the isolated Prague cluster in 2011 and the contraction by 2017, and the temporal check still passes. Its bands for the quarterly totals are wider than rung 1's: the negative binomial admits that a quarter's count is noisier than Poisson.

#9 · Model comparison with LOO

Every rung is a large improvement on the last, by many multiples of dse. The biggest step is adding space at all, the next is letting the map move, and the negative binomial is worth several hundred more - while reducing the effective number of parameters p by almost two thirds. That is the regularisation we saw in the hyperparameters, showing up as predictive accuracy: the Poisson field was fitting bursts as if they were signal.

Now the small print. ArviZ warns about Pareto-$k$ values above 0.7, and the diag_elpd column shows that for the two interaction rungs this is not a handful of points but hundreds of the 3888. That is the price of a very flexible latent field with a short lengthscale: each non-empty cell-quarter has a strong pull on its own neighbourhood, so the posterior without that observation would differ noticeably - the case importance sampling handles badly - and the elpd of rungs 3 and 4 is probably optimistic. Gaps of hundreds and thousands leave the ranking in little doubt, but the numbers should not be quoted to the last digit. The remedy is K-fold cross-validation (az.loo_kfold), ideally leaving out whole blocks of neighbouring cells; at a minute per fit we skip it here, because we have something better than any in-sample estimate: data the models have never seen.

#10 · Forecasting the held-out years

No new machinery is needed. The innovations of 2018-19 - eight steps of $g$ and two yearly increments of the field - were sampled along with everything else, from their priors given each draw's sigma_g, eta_d and $\ell$. simulate already returned all 44 quarters; so far we only looked at the first 36.

Here the total/share parametrisation pays off a second time. In the textbook form $\alpha + g_t + f_t(s)$, a field whose variance grows with the forecast horizon inflates the expected total by $e^{\text{var}/2}$: the forecast number of earthquakes would rise purely because the model is unsure where they will be. With normalised shares the total is governed by $g_t$ alone.

Read the table and the figures together.

What the forecast gets right.

• Where. The hold-out log score improves at every step from rung 1 to rung 3, in both years. The separable model forecasts with the nine-year average map, which still lights up areas that had gone quiet; the interaction model forecasts with the latest map, blurred by one and then two years of random-walk increments, and most held-out epicentres (blue) fall on it. LOO's ordering of rungs 1-3 is confirmed by genuinely unseen data.
• Calibration of the totals. Seven of the eight held-out quarterly totals lie inside rung 4's 90% band, and the eighth (the 7 events of 2019 Q2) is just below it.

What deserves less credit.

• Several small clusters of held-out events lie in the south-west of the window, on cells that both maps leave blank. A random walk can blur the last map; it cannot anticipate a new hot-spot.
• The 2019 map of rung 4 is paler than its 2018 map although the forecast total is the same. That is not a forecast of fewer earthquakes. The panels show the median of each cell, and as the model grows less sure where events will be, each cell's distribution spreads out and its median sinks. Summaries of skewed forecasts need care.
• Cell-level coverage of 99% for a nominal 90% is not a triumph: most cells are empty and an interval that starts at zero covers them trivially. The non-empty cells are the honest number, and there rungs 1, 3 and 4 are about right to slightly wide, while the separable model is too narrow - confident about the wrong map.
• Rungs 3 and 4 are practically tied on the hold-out, although LOO separated them by hundreds. The negative binomial earned its keep in the peak years, full of aftershock bursts; in a quiet period with a few events per cell there is little overdispersion to capture. LOO answers "which model describes data like the training data?", and 2018-19 is not like the training data.
• The band for the total is wide - by late 2019 it spans two orders of magnitude - and the median forecast is flat at roughly the 2017 level, while the truth kept falling: the observed two-year total (208) is less than half of every model's median forecast, and the last quarters sit at the bottom of the band. A random walk has no notion of momentum: "it will stay where it is, give or take a lot" is the only forecast it can make.

#What this model cannot know

If the decline after 2015 was caused by the limits on injection, its cause lies outside the catalogue. The model has no covariate for injection: it sees a rate that went up and then down and extrapolates ignorance in both directions. Had injection resumed in 2018, the forecast would have been exactly the same. It is a descriptive, phenomenological model - excellent for mapping where and when the rate was elevated, with honest uncertainty, and for short-range "more of the same" forecasts of where the next events will be. It is not a tool for deciding what happens if policy changes. For that, the next step is mechanistic: put injection volumes (lagged, and perhaps diffused in space) into $\log\lambda$ as a covariate, and let the GP absorb only what the physics leaves unexplained.

#Try it yourself

1. Magnitudes as marks. The Gutenberg-Richter law says magnitudes above a completeness threshold are exponentially distributed: $m - 3 \sim \text{Exponential}(b \ln 10)$. Estimate the $b$-value (mind the rounding of magnitudes), then let it vary by year (hierarchically). Did the size distribution change as the rate rose and fell? Combine it with the rung-4 forecast to get the probability of at least one M4.5+ event in the window in 2018.
2. Aftershocks, properly. Replace the negative binomial's anonymous overdispersion by a self-exciting term: add $\kappa \cdot y_{c,t-1}$ (last quarter's count in the same cell, and perhaps its neighbours) to $\mu_{ct}$ - a discrete-time Hawkes/ETAS model. Does inv_phi shrink towards zero? Does the field get smoother still? Quarters are coarse for aftershock sequences: try months.
3. Change the dynamics. (a) Make the field increments mean-reverting (an AR(1) on the coefficients instead of a random walk) and compare hold-out scores. (b) Replace the random walk for $g$ with a local linear trend (a random walk on the slope) so that the forecast can carry the post-2015 decline forward - then ask whether you would have trusted the same model in 2013, when the slope pointed up. (c) If you can obtain injection volumes, add them as a covariate.

Next: the C-series challenges, where nobody tells you which rung to stop at.