{'fraction_non_alphanumeric': 0.04172745841844178, 'fraction_numerical': 0.017508024511234316, 'mean_word_length': 4.237249283667622, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Exploring seismic hazard in the Groningen gas field using adaptive kernel smoothing and inhomogeneous summary statistics

The discovery of gas in Groningen in 1959 has been a massive boon to the Dutch economy. From the 1990s onwards though, gas production has led to induced seismicity. In this paper, we carry out a comprehensive exploratory analysis of the spatio-temporal earthquake catalogue. We develop a non-parametric adaptive kernel smoothing technique to estimate the spatio-temporal hazard map and to interpolate monthly well-based gas production statistics. Second- and higher-order inhomogeneous summary statistics are used to show that the state of the art rate-and-state models for the prediction of seismic hazard fail to capture inter-event interaction in the earthquake catalogue. Based on these findings, we suggest a modified rate-and-state model that also takes into account changes in gas production volumes and uncertainty in the pore pressure field.


Introduction
In 1959, a large gas field was discovered in Groningen, a province in the North of The Netherlands. Its recoverable gas volume has been estimated at around 2, 900 billion Normal cubic metres (Nbcm). The extraction rate has varied considerably over the years. After a modest start, large amounts were being produced annually during the early 1970s rising to about 85 Nbcm in 1976. During the next decade, the production volumes tended to decrease, followed by somewhat higher values during the first half of the 1990s. Production fell again during the second half of the decade, before rising in the new millennium to over 53 bcm in 2013. From the 1990s earthquakes were being registered in the previously tectonically inactive Groningen region. Especially the one near Huizinge in August 2012 with a magnitude of 3.6 led to a public demand for a reduction of gas production volumes. The government reacted with legislation to phase out gas extraction and, by 2020, production had fallen to less than 8 Nbcm.
Numerous studies on the Groningen field have been conducted. For example, Geerdink [14] models the times in between earthquakes in terms of the cumulative and annual production rates, pressure, subsidence and fault zones. More recent examples of such a study include Post et al. [23] and Trampert et al. [26]. Van Hove et al. [16] propose a Poisson auto-regression model for the annual hazard maps in terms of subsidence, fault lines and gas extraction in previous years. Both Hettema et al. [15] and Vlek [27] explore the temporal development of seismicity in Groningen by proposing a linear model for the relation between the number of earthquakes over specific periods and gas production volumes. Sijacic et al. [25] focus on the detection of changes in the rate of a temporal Poisson point process by Bayesian and frequentist methods. Bourne et al. [3] modify Ogata's space-time model [22] to include changes in stress level and estimate the probability of fault failures. Other papers, notably Candela et al. [4], Dempsey and Suckale [10] and Richter et al. [24], discuss the modelling of seismicity in relation to stress changes based on a differential equation and embed these in a space-time Poisson point process.
In a previous paper [2], we investigated the temporal development of seismicity in Groningen including data up to 2020 using cumulative and recent gas production as dependent variables in a regression model. We concluded that a decrease in production leads to decreased seismicity. Here, we extend the analysis to take into account spatial variations. First, we calculate non-parametric estimates for the spatio-temporal hazard map by means of an adaptive kernel smoother (cf. Abramson [1], Davies et al. [9] and Van Lieshout [20]) and generalise the bandwidth selection approach suggested by Van Lieshout [20] to the space-time domain. Using this map, we investigate whether a Poisson point process model would suffice. Employing inhomogeneous summary statistics, we find that there is interaction which cannot be accounted for by Poisson models, including the state of the art rate-and-state models in Candela et al. [4], Dempsey and Suckale [10] and Richter et al. [24]. Since rate-and-state models rely on differential equations for changes in Coulomb stress or, equivalently, pore pressure, we shift our attention to pressure and production data in the public domain. The production values are measured monthly at wells whereas pressure is gauged at irregular times at some wells as well as at several observation and seismic monitoring stations. To obtain a gas production map, adaptive kernel smoothing applies. Since mass must be preserved, Van Lieshout's local edge correction [18] is required. For the pressure values, we fit a Gaussian random field, the mean function of which is modelled as a polynomial in space and time. Finally, we propose a modification of the rate-and-state models of Candela et al. [4], Dempsey and Suckale [10] and Richter et al. [24] that exhibits clustering, accounts for the uncertainty in pore pressure, takes into account the varying gas production, and is amenable to monitoring by means of Markov chain Monte Carlo methods based on the history of recorded earthquakes.
The plan of this paper is as follows. In Section 2, we describe the data. Section 3 carries out a comprehensive second-order analysis, Section 4 is devoted to extrapolation of gas production and pore pressure measurements from wells to field. The paper closes with our proposed modification of the Coulomb rate-and-state seismicity model.

Data
Data on the Groningen gas field and the induced earthquakes is available at various sources.  Figure 1: Spatial (left-most panel) and temporal (right-most panel) projections of the 332 earthquakes of magnitude 1.5 or larger with epicentre in the Groningen gas field that occurred in the period from January 1st, 1995, up to December 31st, 2021.

Shapefiles for the Groningen gas field
Shapefiles for the Groningen gas field can be downloaded from the Geological Survey of the Netherlands TNO website www.nlog.nl/bestanden-interactieve-kaart. The files are updated monthly. In this paper we use the map that was published in April 2022. The coordinates of the field are given in the UTM system using zone 31 with metre as the spatial unit, which we rescale to kilometre. The boundary is outlined in the left-most panel of

Earthquake catalogue
An earthquake catalogue for The Netherlands is being maintained by the Royal Dutch Meteorological Office (KNMI) at www.knmi.nl/kennis-en-datacentrum/dataset/aardbevingscatalogus. Data on the period before 1995 is not reliable due to the inaccuracy of the equipment used. Moreover, a threshold on the magnitude is necessary to guarantee data quality. According to Dost et al. [12], for data from 1995, earthquakes with magnitude 1.5 or larger can be reliably recorded; a threshold of 1.3 can be used for the period from 2010 onwards due to an extension of the monitoring network (cf. Hettema et al. [15]). We use data over the time window 1995-2021 and therefore work with a magnitude 1.5 threshold. The coordinates of the epicentres are listed in terms of latitude and longitude. To avoid distortions and for compatibility with the gas field map, we project them to UTM (zone 31) coordinates. This procedure results in 332 earthquakes, the spatial and temporal projections of which are shown in Figure 1.

Wells
The exploration and production company NAM maintains a number of production, injection and observation wells. Their coordinates (in the Amersfoort projected coordinate system used in The Netherlands) are available from the production plans published on www.nam.nl/gas-en-olie/groningen-gasveld/winningsplan-groningen-gasveld.html. For compatibility, we transform the Amersfoort coordinates to UTM (zone 31).
Of the 52 locations in the Groningen gas field shown in Figure 2, 29 are production wells (indicated by a cross), one is an injection well (indicated by a circle) and 22 are observation wells (indicated by a triangle).
A few remarks are in order. Firstly, some production wells (at Midwolda, Noordbroek, Nieuw Scheemda and Uiterburen) were taken out of production around the year 2010 and are no longer in use. Secondly, two south-westerly observation wells (at Kolham and Harkstede) were drilled in a peripheral field rather than in the main reservoir. Finally, up to the mid 1970s, small amounts of gas were extracted from wells not earmarked for production.

Gas extraction
Monthly production values from the start of preliminary exploration in February 1956 up to and including December 2021 were kindly provided by Mr Rob van Eijs from Shell for all of the 29 production wells. The figures were given in cubic metres which we rescale to Nbcm. In Figure 3, we show the time series for six wells chosen to show a range of production patterns: Bierum in the North-East, De Paauwen in the centre and Eemskanaal in the West of the gas field, Amsweer in the central East, Tusschenklappen in the South-West and Zuiderpolder in the South-East.
One may observe that not all wells were drilled at the same time and that some were not in use during the entire period. Also there are differences in the amount of gas extracted: the production figures for Eemskanaal-13 are lower than average. The sharp decline in production from 2014 following legislation is readily apparent.

Pore pressure observations
On nam-feitenencijfers.data-app.html/gasdruk.html, pore pressure observations are available over the period from April 1960 until November 2018. In total, there are 2056 observations. However, these data need some cleaning, as discussed in the following paragraphs. Missing coordinates Since the coordinates of one of the stations are not listed in the production plans (cf. Section 2.3) and therefore unknown, we omit the corresponding 28 pore pressure measurements from consideration. We also disregard the two measurements from an observation well located outside the Groningen gas field. The Eemskanaal-13 well is the only one depleting a peripheral field, the so-called Harkstede block. Moreover, as can be seen from Figure 3, it is extracting less gas than other wells. The combined effect is that the pore pressure measurements are somewhat higher than at other wells. According to an expert, setting its location to either that of the Eemskanaal plant or to that of the installation at Harkstede would lead to biases and it is therefore preferable to ignore the seven observations for Eemskanaal-13 altogether.
Invalid measurements The NAM file mentions ten cases in which the observations are invalid for various reasons. We delete these measurements.
After cleaning, we are left with 2009 pore pressure measurements. Figure 4 shows time series at four locations: a production well in the South-West (Slochteren), two observation wells (Harkstede and Kolham) in peripheral fields and the injection well at Borgweer in the East. Mostly, the graph is initially flat, followed by a decrease. Note that the measurements in the periphery are higher than those in the main reservoir.

Exploratory data analysis
In this paper, we treat the earthquake catalogue as a spatio-temporal point pattern, a realisation of a point process in space and time. Formally, let Ψ be a simple spatio-temporal point process in W S × W T for bounded open sets W S ⊂ R d and W T ⊂ R and suppose that its first order moment measure exists, is finite and absolutely continuous with respect to the product of Lebesgue measures ℓ in space and time (see e.g. Chiu et al. [6]). Write λ for its Radon-Nikodym derivative, known as the intensity function, and N (Ψ ∩ (W S × W T )) for the number of points placed by Ψ in W S × W T . Intuitively speaking, λ(s, t)dsdt is the probability that Ψ places a point in the infinitesimal region dsdt around (s, t) ∈ W S × W T .
As usual in spatial statistics, we start by an empirical exploration of the trend and interaction.

Adaptive kernel estimation of the intensity function
Our first step is to estimate the intensity function of the point process of earthquakes. The standard technique to do so is kernel estimation as proposed by Diggle [11]. This technique can be seen as a generalisation of the histogram. Briefly, a Gaussian kernel (say) with a given standard deviation h S in space and h T in time is centred at each of the points in a realisation of Ψ and their sum is reported. The choice of the bandwidhts h S and h T is crucial. Some asymptotical results are available in Chacón and Duong [5], Van Lieshout [19] and Lo [21], but practical rules of thumb seem to be lacking.
The risk of using the same bandwidths h S , h T at all points of Ψ is that, as all onesize-fits-all solutions, this approach tends to oversmooth in regions that are rich in points while at the same time it does not smooth enough in sparser regions. To overcome this drawback we propose to use an adaptive smoother as introduced for classic random variables by Abramson [1]. In a spatial context, such estimators were studied by Davies et al. [9] for Poisson point processes and in Van Lieshout [20] for general point processes with interaction between the points. The underlying idea is to weigh the bandwidth at (s, t) ∈ Ψ by a scalar that is inversely proportional to λ(s, t). The power 1/2 is motivated by asymptotics (cf. Abramson [1] and Van Lieshout [20]), in practice other powers could also be used. Formally, set Here H(h 1 , h 2 ), h 1 , h 2 > 0, is a 3 × 3 diagonal matrix whose first two entries are h 1 and whose third entry is h 2 , w((s, t), h S , h T ) is an edge correction term and κ is a symmetric probability density function, the kernel. Since λ is unknown, c(s, t) must be estimated. We will do so by plugging in a standard kernel estimator with fixed bandwidth.
To select the bandwidth, we propose the following completely non-parametric and computationally easy algorithm.
Algorithm 1. Let ψ be a non-empty spatio-temporal point pattern that is observed in the window W S × W T . Then 1. Choose global bandwidth h g,S , h g,T by minimising with local edge correction weights 3. Choose adaptive bandwidths h a,S , h a,T , by minimising over h S , h T > 0 (with minimal h 2 S h T in case of multiple solutions) whereλ A is given by (1) and (2) with w ≡ 1 upon plugging in the edge-corrected pilot estimatorλ c for λ.
Local edge correction weights, as suggested by Van Lieshout [18], for the adaptive kernel estimator take the form h 2 a,S h a,T (ĉ(y)) 3 dz and ensure that (1) is mass preserving in the sense that its integral over W S × W T is equal to the number of points in ψ ∩ (W S × W T ). In selecting the bandwidth in steps 1 and 3 of Algorithm 1, no edge correction is applied in order to obtain a clear optimum (see Cronie and Van Lieshout [8]). The justification of Algorithm 1 lies in the Campbell-Mecke theorem (cf. Chiu et al. [6]), which states that The first step is a modification to space-time of the Cronie-Van Lieshout bandwidth selector [8] for purely spatial point processes. An important difference is that in space one optimises over a single parameter and usually, but not always, there is only one minimiser. In the space-time domain, the optimisation is done with respect to two parameters, h S and h T . Therefore, as a rule, a curve of minimisers is found in steps 1 and 3 of Algorithm 1. We pick the optimiser having the smallest scale h 2 S h T . For our earthquake catalogue, the results are shown in Figures 5 and 6. The pilot bandwidths are h g,S = 9.4 and h g,T = 182.5. The corresponding edge-corrected projections on space and time are shown in, respectively, the left-most panel of Figure 5 and the broken line in Figure 6. They should be compared to the projections of the edge-corrected adaptive kernel estimate with h a,S = 6.9 and h a,T = 212.9 shown in the right-most panel in Figure 5 and the solid line in Figure 6. Note thatλ A attains higher values in the central reservoir thanλ, lower values near the eastern border of the gas field. From a temporal perspective, the extremes in years with a large number of earthquakes are somewhat more pronounced for the adaptive kernel estimate.

Inhomogeneous spatio-temporal K-and J-functions
Having estimated the trend, we now turn our attention to the inter-point interactions. To quantify such interaction, information about the joint distributions of pairs, or tuples, of points is required, which is formalised by the higher-order analogues of the intensity function, the product densities λ (n) . Heuristically, λ (n) ((s 1 , t 1 ), . . . , (s n , t n ))ds 1 dt 1 · · · ds n dt n is the probability that Ψ places points at each of the infinitesimal regions ds i dt i around (s i , t i ), i = 1, . . . , n.
A spatio-temporal point process Ψ on R 2 × R is said to be intensity-reweighted moment stationary (IRMS) (cf. Cronie and Van Lieshout [7]) if its product densities λ (n) of all orders exist,λ = inf (s,t) λ(s, t) > 0 and, for all n ≥ 1, ξ n is translation invariant in the sense that ξ n ((s 1 , t 1 ) + (a, b), . . . , (s n , t n ) + (a, b)) = ξ n ((s 1 , t 1 ), . . . , (s n , t n )) for almost all (s 1 , t 1 ), . . . , (s n , t n ) ∈ R 2 × R and all (a, b) ∈ R 2 × R. Here ξ n are the n-point correlation functions defined in terms of the λ (n) by setting ξ 1 ≡ 1 and for other n recursively by where D 1 ,...,D k is a sum over all possible k-sized partitions {D 1 , . . . , D k }, D j = ∅, of the set {1, . . . , n} and |D j | denotes the cardinality of D j . For a Poisson point process, in which there are no correlations between the points, ξ n ≡ 0 for n ≥ 1. The point process is said to be second order intensity-reweighted stationary (SOIRS) if the translation invariance holds up to n = 2 (see Gabriel and Diggle [13]). Various summary statistics exist to explore inter-point interactions. Write |t| ≤ h T } and let Ψ be an IRMS spatio-temporal point process. Then the statistic t 1 ), . . . , (s n , t n )) n i=1 ds i dt i quantifies cumulative correlations of order n = 1, 2, . . . up to ranges h S ≥ 0 in space and h T ≥ 0 in time. Multiple orders can be combined. For example, Van Lieshout [17] proposed for all spatial ranges h S ≥ 0 and temporal ranges h T ≥ 0 for which the series is absolutely convergent. Values greater than one are indicative of repulsion, values smaller than one suggest clustering. For further details, the reader is referred to Cronie and Van Lieshout [7]. Truncating at n = 2, in terms of the inhomogeneous K-function [ξ 2 ((0, 0), (s 1 , t 1 )) + 1] ds 1 dt 1 of Gabriel and Diggle [13], which is well-defined under the SOIRS assumption. Values greater or smaller than the volume of S h T h S are indicative of, respectively, clustering and inhibition between points. In practice, one must estimate these statistics based on a pattern observed in W S × W T . The definition of the inhomogeneous J-function does not immediately suggest a suitable estimator. However, Cronie and Van Lieshout [7] showed that, for all y ∈ R 2 × R, whenever well-defined and the denominator is greater than zero. Here the function u y h S ,h T is defined by and G is the generating functional of Ψ, that is, for h S , h T ≥ 0, under the convention that empty products take the value one. G !y is defined similarly in terms of the distribution of Ψ\{y} given there is a point of Ψ at y. Being expectations, the numerator and denominator in equation (6) can be estimated in a straightforward manner. Write W ⊖h S S for the set of points in W S that are at least h S away from the border of W S and let W ⊖h T T be the similarly eroded temporal domain. Then, given a finite point is an unbiased estimator of the denominator in (6). An unbiased estimator for the numerator is obtained analogously. Finally, is an unbiased estimator of K inhom (h S , h T ). Note that the intensity functions are unknown. A practical solution is to plug in their estimated counterparts (cf. Section 3.1).
For the earthquake catalogue depicted in Figure 1, consider Figure 7. The solid line in the left-most panel of the figure is the graph of (K inhom (r, 100r)/(2π)) 1/3 for the earthquake data. It lies above the graph of the same function for a Poisson point process (shown as a broken line in the left-most panel of the figure), suggesting attraction between the points. In the right-most panel of Figure 7, the graph ofĜ(1 − u y r,100r ) estimated from the data lies above that of G !y (1 − u y r,100r ), which confirms the suggested clustering.

Explanatory variables
Gas production and pore pressure (cf. Section 2.4 and 2.5) may well have an effect on the earthquake rates. However, they are measured only for a limited number of locations and times. Therefore, in this session, we discuss how to calculate appropriate values for the entire space-time domain.

Non-parametric smoothing of gas production values
Monthly production values are available for the 29 production wells shown by a cross in Figure 2 over the time period 1956-2021. The gas extracted from the Eemskanaal-13 bore hole is listed separately from the other ones at the Eemskanaal site. The reason is that the Eemskanaal-13 pipe is not drilled vertically but bends underground in such a way that it is depleting the Harkstede block (cf. Section 2.5). Therefore we assign its production to the location of the observation well at Harkstede. Thus, we end up with a spatial pattern of 30 points, the locations of the 29 production wells as well as the Harkstede proxy for the Eemskanaal-13 pipe. In order to smooth the production out in a mass-preserving fashion, an adaptive kernel method similar to that proposed in Section 3 for the intensity function can be used. More specifically, write V (s, t) ≥ 0 for the volume of gas produced at a well at s in month t and write φ for the well pattern. We then use the spatial counterpart of Algorithm 1 proposed by Van Lieshout [20] to select bandwidths h g,S and h a,S and set, for s 0 ∈ W S , t 0 ∈ W T , based on spatial kernel κ and using local edge correction weights w, cf. Section 3. In time, V (s, t) is spread evenly over the ℓ(m(t)) days in the month m(t) to which t ∈ W T belongs. This smoother is mass preserving, as For the data discussed in Section 2.4, Algorithm 1 yields a pilot bandwidth h g,S = 7.6 and h a,S = 6.9. The resulting gas production maps for January 2012, before legislation to phase out gas extraction came into effect, and for January 2021 are shown in Figure 8. The decrease in extracted volume is evident. Moreover, the remaining production is located mostly in the southern part of the gas field.

Regression analysis for pore pressure
From Figure 4, it is clear that the pore pressure in the reservoir is decreasing in time. The trend is similar for most wells, with some exceptions, for example because of large faults or at the periphery of the field. Recall from Section 2 that our earthquake catalogue contains tremors from January 1st, 1995 onward. If we wish to use the pore pressure as an explanatory variable in a monitoring model, we need their values over the same time period. Thus, our goal in this section is to perform a regression analysis using only the 352 pore pressure measurements from January 1st, 1995, and later. Suppose that the observed pore pressures can be seen as realisations of stochastic variables X(s, t) that can be decomposed in a trend term m((s, t); β) and noise E(s, t) as follows: X(s, t) = m((s, t); β) + E(s, t).
We assume that the E(s, t) are independent and normally distributed with mean zero and variance σ 2 . For the trend fit a polynomial in space and time, where we centre the spatial locations at s 0 = (750, 5900) in the UTM system (zone 31 in units of km) and take days as the temporal unit counting from January 1st, 1995. The parameter β = (β 1 , . . . , β 26 ) can be estimated by the least squares method. The fitted pore pressure maps at midnight January 1st in the years 1995 and 2022 are given in Figure 9. One may note the elevated values in the Harkstede block in the South-West as well as those in the Northern border regions. To validate the model, the actual pore pressure measurements are plotted against their fitted values in the left-most panel of Figure 10. The graph seems reasonably close to a straight line. The histogram of the residuals is shown in the right-most panel of Figure 10.
Most residuals are quite small and the histogram is centred around zero with estimated standard deviationσ = 7.17 bara.

Discussion and further work
In this paper, we carried out an exhaustive second order exploratory analysis of the spatiotemporal point pattern of earthquakes recorded in the Groningen gas field since January 1995. To do so, we needed to develop new methodology. We proposed an adaptive kernel smoothing technique for estimating the intensity function and suggested a practical algorithm for selecting the spatial and temporal bandwidths. The estimated intensity function was then plugged into state of the art inhomogeneous summary statistics to quantify the degree of clustering in the earthquake catalogue. We also applied our new adaptive kernel smoothing technique to monthly gas production figures. Finally, we performed a regression analysis on pore pressure data for the gas field.
In the rate-and-state models (e.g. Candela et al. [4], Dempsey and Suckale [10] and Richter et al. [24]) that can be seen as the state of the art in modelling the seismic hazard and that is being used for planning, the earthquake intensity λ (the rate) is assumed to be inversely proportional to a state variable Γ, that is, The state variable Γ(s, t) is defined by the ordinary differential equation where S is the Coulomb stress with a reduced friction coefficient and α > 0. According to Richter et al. [24], Coulomb stress changes are proportional to changes in pore pressure X(s, t), that is dS(s, t) = −βdX(s, t) for some scalar β > 0. Consequently, the differential equation for the state can be written as dΓ(s, t) = α [dt + βΓ(s, t)dX(s, t)] .
In the rate-and-state model, the earthquakes constitute a Poisson point process with intensity function λ, possibly modified by a fault map. The parameters α and β and the initial state Γ(s, 0) ≡ 1/λ 0 are treated as unknowns and can be estimated, for example by the maximum likelihood method.
Based on the exploratory analysis in Sections 3 and 4, the rate-and-state model can be criticised on several points. As we saw in Section 3, the earthquake pattern exhibits clustering. Since by definition the points in any Poisson point process do not interact with one another, the apparent clustering of earthquakes cannot be described by the rate-and-state model. Secondly, the pressure values are assumed to be known everywhere, in practice by interpolation of the measurements. Proceeding in this way, the uncertainty in the interpolations is ignored. It would be better to treat the X(s, t) as a random field. Lastly, the varying gas extraction is not taken into account.

Γ(s, t)
whereṼ (s, t) is the gas extracted at s during the year preceeding time t and U (s, t) is a correlated Gaussian field that accounts for random effects. Monitoring can then be based on the posterior distribution of Λ or, equivalently Γ and U , given the recorded earthquakes. The implementation requires careful use of Markov chain Monte Carlo techniques and is the topic of our ongoing research.