Statistics of seismic events at the Groningen field

Abstract

Depletion of gas fields, even in a tectonically inactive area can induce earthquakes. This is the case for the Groningen gas field, located in the north of the Netherlands. Increased seismic activity raised public concern which led to the government trying to understand the cause of the earthquakes and optimize production such as to minimize the risk of induced seismicity. The main question is how production is correlated with induced seismicity. In this paper we deal with the statistics of seismic events using Bayesian model comparison and a Bayesian change point model. We have developed a method to assess seismic event rate, its changes and tendencies. These statistical analyses are in agreement with each other and find a constant event rate up to 2003, an increasing event rate from 2003 to 2014 and a preference for a decreasing event rate from early 2014 to now. Seasonality in the production and the number of events is present. The seasonality indicates a delay ranging between 2 and 8 months between seismicity and production changes. The question of interest is whether the production reduction since January 2014 has had an effect on the seismicity occurring in the Groningen field. The number of events in the Groningen field in the area affected by the production change has been reduced significantly. We present evidence that changes in seismicity are closely related to changes in production.

Appendix: Normalization of cross correlations

The cross-correlation of two functions f and g on the real t-axis is defined as the Riemann–Stieltjes integral

$${\text{h}}\left( {\text{k}} \right) = \int {{\text{f}}\left( {\text{t}} \right){\text{g}}\left( {{\text{t}} + {\text{k}}} \right)\;{\text{d}}\alpha \left( {\text{t}} \right)}$$

Either we choose α = t or take a step function for α. In this way we incorporate (Riemann) integrals and series into one formalism (Rudin 1976, Chapter 6). We leave the integration boundaries unspecified as yet.

Now form ${\text{Q}} \equiv \int {\left[ {{\text{f}}\left( {\text{t}} \right) - \lambda {\text{g}}\left( {{\text{t}} + {\text{k}}} \right)} \right]^{2} {\text{d}}\alpha \left( {\text{t}} \right)}$, λ being a real number. Then Q is obviously non-negative.

Hence $\lambda^{2} \int {{\text{g}}\left( {{\text{t}} + {\text{k}}} \right)^{2} {\text{d}}\alpha \left( {\text{t}} \right)} - 2\lambda \int {} {\text{f}}\left( {\text{t}} \right){\text{g}}\left( {{\text{t}} + {\text{k}}} \right){\text{d}}\alpha \left( {\text{t}} \right) + \int {{\text{f}}\left( {\text{t}} \right)^{2} {\text{d}}\alpha \left( {\text{t}} \right) \ge 0}$.

In order that this quadratic equation in λ satisfies this inequality for all functions f and g the discriminant should be non-positive. This leads to the general requirement that

$$\left| {\int {{\text{f}}\left( {\text{t}} \right){\text{g}}\left( {{\text{t}} + {\text{k}}} \right){\text{d}}\alpha \left( {\text{t}} \right)} } \right|^{2} \le \int {} {\text{ f}}\left( {\text{t}} \right)^{2} {\text{d}}\alpha \left( {\text{t}} \right)\int {{\text{g}}\left( {{\text{t}} + {\text{k}}} \right)^{2} {\text{d}}\alpha \left( {\text{t}} \right)}$$

If we now set the integration boundaries to −∞ and +∞ we retrieve the familiar Cauchy–Schwarz inequality, since k can then be omitted from the last integral for the given choices of α. We can now normalize the cross-correlation with the square root of the auto-correlations of f and g taken at lag 0, with the neat result that the cross-correlation so normalized is confined to [− 1, 1].

However, we are obviously dealing with data over a finite time interval. We normalized the cross-correlation function with the square root of the right hand side, taking k = 0, since lag-dependent normalization gives obviously rise to distortion, which is unacceptable. But then we have no guarantee that the cross-correlation so normalized always remains between −1 and 1.