Stress-based, statistical modeling of the induced seismicity at the Groningen gas field, The Netherlands

Groningen is the largest onshore gas field under production in Europe. The pressure depletion of the gas field started in 1963. In 1991, the first induced micro-earthquakes have been located at reservoir level with increasing rates in the following decades. Most of these events are of magnitude less than 2.0 and cannot be felt. However, maximum observed magnitudes continuously increased over the years until the largest, significant event with ML=3.6\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_L=3.6$$\end{document} was recorded in 2014, which finally led to the decision to reduce the production. This causal sequence displays the crucial role of understanding and modeling the relation between production and induced seismicity for economic planing and hazard assessment. Here we test whether the induced seismicity related to gas exploration can be modeled by the statistical response of fault networks with rate-and-state-dependent frictional behavior. We use the long and complete local seismic catalog and additionally detailed information on production-induced changes at the reservoir level to test different seismicity models. Both the changes of the fluid pressure and of the reservoir compaction are tested as input to approximate the Coulomb stress changes. We find that the rate-and-state model with a constant tectonic background seismicity rate can reproduce the observed long delay of the seismicity onset. In contrast, so-called Coulomb failure models with instantaneous earthquake nucleation need to assume that all faults are initially far from a critical state of stress to explain the delay. Our rate-and-state model based on the fluid pore pressure fits the spatiotemporal pattern of the seismicity best, where the fit further improves by taking the fault density and orientation into account. Despite its simplicity with only three free parameters, the rate-and-state model can reproduce the main statistical features of the observed activity.


Introduction
Induced seismicity has been observed as a consequence of various human underground activities for energy production and resource extraction, such as (1) fluid or gas removal in the subsurface, (2) mining or quarry sites, (3) fluid or gas injection and (4) water impoundment sites (Foulger et al. 2018;Evans et al. 2012;Ge et al. 2009). However, the causal relationship between underground activity and seismicity is complex, as at several places induced seismicity is absent or only of small magnitude, while in other cases, independent of the activity type, large earthquakes have occurred, including examples like the 1976 M = 7.0 Gazli, Usbekistan (Simpson and Leith 1985), 1989 M L = 5.6 Völkershausen, Such large earthquakes increased the awareness of induced seismicity worldwide, especially if the events occurred in regions of low tectonic seismicity (Keranen et al. 2014;van Thienen-Visser and Breunese 2015). As a consequence, ways of regulating this process are increasingly discussed (Langenbruch and Zoback 2016) or imposed (van Thienen-Visser and Breunese 2015). Such regulations are often based on vague assumptions or empirical observations with poor statistics. Therefore, theoretical models should be employed to fill this knowledge gap. These must be developed and verified to reliably forecast induced seismicity as a function of a few input parameters.
According to the current physical models, induced earthquakes occur on pre-existing faults in response to strength and stress changes. In particular, changes of the Coulomb-Failure stress (CFS), the difference between the shear stress and the product of the friction coefficient and compressive normal stress, are known to be a first-order quantification of the criticality of faults (Toda et al. 2005). Although the general mechanism seems to be clear, the triggering of seismicity depends on many parameters, such as the distribution and orientation of pre-existing faults, the state of effective stress resolved on the faults in relation to their frictional strength, as well as the acting stressing rate. All these parameters are difficult to measure. For instance, the existence of blind faults in a low-seismicity region is often only recognized after a major induced or triggered earthquake. As well, the state of stress on faults is often unknown. Differences in the critical stress threshold needed to trigger the first earthquakes indicate large variability for different regions, even if the production technique and intensity are similar (van Eijs et al. 2006).
While, e.g., pressure perturbations of ∼ 0.1 MPa triggered earthquakes in Oklahoma in Central USA, a reduction of the reservoir pore pressure by ∼ 10 MPa was necessary to trigger earthquakes in the Groningen gas field in The Netherlands Fig. 1 a Location map of Groningen and other gas fields in Northern Germany with observed seismicity. b Yearly changes in pressure (black) and compaction (blue) at the Groningen gas field and the seismicity with M ≥ 1.5 recorded by KNMI (Candela et al. 2018). In our study we use the Groningen gas field as test case. Figure 1 shows the temporal evolution of the induced seismicity over the production time.
The clearly delayed onset of the seismicity in Fig. 1b has been explained by a subcritical state of stress on the faults in the reservoir. Among others, Dempsey and Suckale (2017) assumed a linear relation between the pore pressure change and the seismicity rate on isolated faults in Groningen, considering the stress threshold as an additional free parameter to model the long delay of seismicity. Using this method, the estimation of the detailed distribution of the pre-stress levels within the fault network can be only done retrospectively and involves a number of free parameters.
The linear Coulomb failure model combined with a stress threshold has a consistency problem. It implies a zero effective tectonic stressing and thus a completely aseismic environment before the critical stress is reached. However, we know from many examples that intraplate regions also experience continuous tectonic stressing, even if the absolute rates are small. Once the threshold is reached this model results in a non-constant earthquake rate, even without human interference. An alternative model, commonly known as Coulomb rate-and-state (RS) model, was developed by Dieterich (1994), who derived a constitutive law for the evolution of the seismicity rate on a population of faults governed by rate-and-state-dependent friction, subject to a stress perturbation (see "Seismicity model and testing approach"). The RS model has already been widely applied to model seismic activity, e.g. associated with mainshocks (Hainzl et al. 2009;Cattania et al. 2015), magma intrusions (Toda et al. 2002), water dam loading , waste water disposals (Norbeck and Rubinstein 2018), and most recently for subregions of the Groningen gas field (Candela et al. 2019). The RS model assumes a constant tectonic loading leading to a constant low seismicity rate in the absence of any human interference and does not involve a threshold level for triggering. However, the RS model is able to explain the delayed onset of the seismicity, depending on the time function of the anthropogenic stress loading.
Another key question to be answered is what parameter controls the rate of induced seismicity. In a general analysis, van Eijs et al. (2006) found three key parameters for hydrocarbon reservoirs in the Netherlands: a critical pressure drop, a stiffness contrast between seal and reservoir rock, and the fault density. For the Groningen gas field different suggestions have been made. Shapiro (2018) found a linear relationship between produced gas volume and seismicity after the onset of the seismicity and assumed the extracted fluid volume rates as the input in their modeling approach. Dempsey and Suckale (2017) took the pore pressure changes as input data to calculate the Coulomb failure stresses on fault segments based on poroelastic effects in a simplified fault system (see Fig. 3). Similarly, Candela et al. (2019) based their model on the pore pressure combined with detailed Coulomb failure stress calculations on the faults for selected regions in the reservoir and the non-linear rateand-state approach. In contrast, Bourne et al. (2014) based their model on the vertical compaction data to calculate the compaction strain as the driving mechanism for seismicity. This model was modified in their latest publication (Bourne et al. 2018), where the authors include the reservoir pore pressure and the strain as supposed by van Wees et al. (2018) and added the topographic gradient of the reservoir layer combined with an extreme threshold failure model to predict the seismicity. Buijze et al. (2017) indicated the fault offset as another controlling parameter.
In this paper, we explore the RS and the Coulomb failure threshold model for explaining the spatiotemporal evolution of the field seismicity observed in Groningen, since the start of production, i.e. including the delayed onset and the change of seismicity during recent years. We analyze and compare different input parameters in statistical tests to evaluate whether pore pressure reduction, reservoir layer compaction, layer thickness or fault distribution and orientation have a sufficient sensitivity to justify their usage in seismicity models.

Coulomb rate-and-state (RS) seismicity model
The most common assumption is that static stress changes lead to earthquake nucleations. Seismicity is promoted when the Coulomb stress increases and inhibited when it decreases. The Coulomb stress is defined as Here is the shear stress in the direction of slip, f is the coefficient of friction, and eff is the effective compressive normal stress on the fault plane; the latter is given by eff = n − p , where n is the Cauchy stress acting on the rock skeleton and p the change in pore pressure in the pore volume of a porous rock.
Based on constitutive friction laws derived in laboratory experiments, Dieterich (1994) derived a constitutive law for the evolution of seismicity on a population of faults governed by rate-and-state-dependent friction, subject to a stress perturbation. The seismicity rate R is given as a function of the stressing history where r is the background rate of seismicity, ̇ is the background shear stressing rate and is the state variable that (2) R = ṙ , 252 Page 4 of 15 depends on the stress evolution. If normal stress changes are small compared to the absolute values, the evolution of is given by Here A and are positive non-dimensional constitutive parameters (Linker and Dieterich 1992). Thus, the input parameter S in the RS model is simply the Coulomb stress with a reduced friction coefficient. However, the earthquake nucleation in a rate-and-state framework is based on physically different principles to the simple Coulomb failure criterion, which is employed in a linear Coulomb model (CM). The RS model employs a frictional instability leading to rupture nucleation that considers the dependence of friction on slip and slip rate. While the CM assumes instantaneous triggering if a certain threshold of the Coulomb stress is exceeded, the nucleation process in the RS model leads to a delayed response in the seismicity. Nonetheless, the concepts underlying Coulomb stress transfer can still be used to interpret the spatial distribution of seismicity over a longer period, because the total number of events triggered in the RS model by a transient stress change is known to be proportional to CFS = CFS(t + t) − CFS(t) (Dieterich 1994;Hainzl et al. 2010a).
In the Groningen gas field, the dominant rupture mechanism is normal faulting (Willacy et al. 2018), indicating that the largest principle compressive stress ( 1 ) is vertical and the smallest ( 3 ) is horizontal. From previous works, it is known that most pre-existing faults dip at around 70 • (van Wees et al. 2018). The gas extraction in Groningen resulted in pressure drops, average around 22 MPa, in the reservoir layer. As a consequence, the horizontal stress h decreases according to poroelasticity with declining reservoir pressure ( p < 0 ) in a linear way, namely with Biot-Willis coefficient and the drained Poisson's ratio (Segall and Fitzgerald 1998;Dempsey and Suckale 2017). This leads to Coulomb stress changes which are proportional to the pore pressure changes In an alternative approach, the reservoir is approximated as a thin sheet, meaning that its lateral extents are significantly larger than its vertical extent. In this approximation, the deformations induced by pore fluid pressure, as an internal body force, will be approximately uniaxial due to symmetry. That means that all strain components related to that change are approximately zero inside the reservoir with the exception of the vertical strain zz which can be approximated by with h being the thickness of the layer (Bourne et al. 2014). This estimation leads to a simple proportionality between the stress changes and compaction where c 2 depends on the dip angle and the assumed seismic coupling coefficient and K is the bulk modulus.
In the following, we test both alternative models. In particular, we assume an induced stressing rate where dots refer to time derivatives.

Input data
Groningen gas field is located in a region of very low tectonic seismicity. No seismicity was recorded in Northern Netherlands until December 1986 when near Assen an earthquake of M L = 2.8 occurred. This event took place nearby to a gas field in production, with a shallow depth (about 1 km). Thus it was classified as an induced event (van Eck et al. 2006). Since 1995 a seismic network of 19 permanent stations (11 borehole stations in the North of the Netherlands) and 17 additional accelerometers has been installed in the Netherlands and recorded several hundred induced events in the northern Netherlands (points in Figs. 1 and 2). The prominent characteristics of the observed seismicity at Groningen gas field are as follows: the first recorded earthquake was a M L = 2.4 event on December 5, 1991, at a depth of 3 km. It happened nearly 30 years after pressure depletion started. The number of events increased in the following years almost exponentially. Similarly, the observed maximum magnitude of the earthquakes continuously increased, until the so-far largest earthquake with M L = 3.6 happened on August 16, 2012. In 2014, the production was drastically reduced and limited, immediately followed by a reduction in seismicity, visible as less dense points in Fig. 1 and shown by the yearly sum in Fig. 5.
We used the seismic catalog from KNMI (Koninklijk Nederlands Meteorologisch Instituut) as it is provided on the KNMI website (www.knmi.nl, events after June 2017, downloaded on January 9, 2019) and as it was used by Dempsey and Suckale (2017) (until June 2017). Since 1995, the magnitude of completeness for the seismic network is 1.5 (Dost et al. 2012); therefore, we restrict the analysis on events with M L ≥ 1.5.
We determine the seismicity as a function of the pore pressure drop in the Groningen reservoir and the (6) zz = h∕h reservoir-layer compaction, both provided by the NAM (Nederlandse Aardolie Maatschappij) for the Groningen gas field on a yearly basis with a spatial resolution of a 500 m × 500 m grid at locations (x n , y n ) with n = 1, … , N and N = 5066 . Production in Groningen was planned to minimize spatial variations in reservoir overpressure. This was possible because the high permeability in the central field ensures the pressure reaches equilibrium within months (van Thienen-Visser and Breunese 2015). However, most production wells are located in the southern part of the reservoir, which may lead to a temporary disequilibrium during enhanced production phases and cause slightly higher pressure drops in the southern area. Nearly all over the gas field the total pressure drop by 2014 ranges from 20 to 28 MPa (Fig. 2).
The reservoir layer compaction model by NAM is based on observed surface subsidence and strain measurements in boreholes. While the pressure drop is almost constant over the field's dimension, the compaction is largest in the western central part of the field (Fig. 2). In this area, the density of faults cutting the reservoir layer is also increased. Also noteworthy is that the reservoir compaction smooths out at the boundaries of the gas reservoir, similar to the displacement solution of a crack model, while the pressure drop model has sharp boundaries including a narrow valley-like structure (reservoir low) bounded by two distinct NW-SE striking faults systems in the western central field. As already recognized by Bourne et al. (2014) and Dempsey and Suckale (2017) most of the induced earthquakes occur in this central western part of the field. For the analysis, we linearly interpolate both pressure and compaction data and calculate their time derivatives in time steps t at each of the spatial grid nodes, For the bulk modulus, we use a value of K = 10 GPa.
To determine the compaction-induced stress, which depends on vertical strain, both the absolute reservoir thickness h(x n , y n ) and its change (compaction) are needed. While the compaction is provided directly from NAM, we estimate h(x n , y n ) from the layer model provided by NAM. The thickness of the Rotliegend formation is taken as the reservoir thickness. The Rotliegend layer gradually thickens towards NW with a mean value of 225 m (Fig. 3). Additionally, the layer dips towards the NW. The reservoir thickness may also influence a depth-normalized seismicity rate. However, due to its small relative changes, we use a constant thickness for both, the integration of seismic production and the stressing of the reservoir layer in model B, Ṡ (t k , x n , y n ) = c 2ḣ (t k , x n , y n )∕h(x n , y n ) (compare Eq. 8).
In our analysis, we additionally use fault density. The fault information provided by NAM are summarized by Dempsey and Suckale (2017). The largest 325 faults are given as vector lines (see Fig. 3). For the fault density map, we set the point spacing describing the faults to 100 m. These points are smoothed using a Gaussian function with a standard deviation of 1 km and the final distribution, (x n , y n ) is normalized. Besides the model with spatially constant background rate r(x n , y n ) = r 0 , we test a model where the background rate is proportional to the fault density, r(x n , y n ) = (x n , y n )r 0 . Note that in both cases, the total background rate for the whole reservoir is the same, namely r tot = Nr 0 . Furthermore, we test the influence of the fault orientation . For that purpose, we assume a normal faulting regime with unknown fault dip

Model implementation
Based on the input fields ṗ(t k , x n , y n ),̇h(t k , x n , y n ), h(x n , y n ) and r(x n , y n ) , we first calculate Ṡ (t k , x n , y n ) for models A and B and then iterate the state variable in each grid point (x n , y n ) according to the algorithm provided by Hainzl et al. (2010b). Due to the absence of historic earthquake data and the large uncertainties in stress estimations, we use the corresponding parameters as fitting parameters. In particular, we use as free model parameter r 0 , A and t a = A ∕̇ . Note that we use for convenience t a instead of ̇ . Furthermore, the factor of proportionality c 1 , respectively c 2 , can be simply absorbed by associating A with A ∕c 1,2 and ̇ with ̇∕c 1,2 in the equations and is not an independent free parameter. Assuming that the model parameters are constant in space, our model consists, therefore, of three free parameters which have to be estimated from the data: t a , A , and r 0 . Starting from constant background rate ( = 1∕̇ , the state parameter is iterated in time steps of fifty days, t = 50 d, in each grid point. This leads to an estimated seismicity rate R(t k , x n , y n ,z) on the average depth level of the reservoir. Our input data have no depth resolution. To account for the total rate in each map location, the rate has to be integrated over the seismogenic width which is here assumed to be equal to the reservoir height, yielding To estimate the free model parameters, we use the maximum likelihood method (Daley and Vere-Jones 2003;Hainzl et al. 2010b). For Z earthquakes occurred at times t i and locations i within a given time interval , the logarithmic likelihood value can be determined by where k(i) and n(i) refer to the index of the time bin and spatial volume in which the ith event fall. However, it is important to consider location uncertainties. For that purpose, we estimate the probability that an earthquake is associated with different spatial grid points according to a Gaussian location error with standard deviation . Note that we normalize all resulting weights w n,i , thus the sum over all grid points is one. Then the spatially averaged rate at the occurrence time is calculated by R(t k(i) ) = ∑ N n=1 w n,i R(t k(i) , x n , y n ) . In the case of = 0 , the weight for the spatial volume in which the ith event fall is one, while all other grid points are given a weight of zero. In this case, Eq. (9) is directly used. In all other cases, Eq. (9) is used after replacing R(t k(i) , x n(i) , y n(i) ) by R(t k(i) ) . In particular, in the case of = ∞ , all spatial grid points get the same weight which is equivalent to only considering the timing of the events, but ignoring the exact location of the events within the reservoir. The r 0 value optimizing ln L can be found analytically for given A and t a (Hainzl et al. 2010b). To find the global optimum, we perform a two-dimensional grid search related to the parameters A and t a . The time interval is set to T 0 = 1959 when the production started until the end of our input data T 1 = 2017 , assuming that no M ≥ 1.5 occurred before 1991.

Alternative models
For comparison, we also test three models with three or fewer free parameters. The first one is the stationary Poisson model, R(t k , x n , y n ) = r 0 , with one free parameter r 0 . This model is the basis of time-independent seismic hazard calculations and represents the random and memoryless occurrence of activity with constant rate. The maximum likelihood method yields for the Poisson model the optimized The second model is the Coulomb model (CM), where the seismicity rate is proportional to the induced stressing rate (if positive), i.e. R(t k , x n , y n ) = r 0 + aṠ(t k , x n , y n ) ⋅ H(Ṡ(t k , x n , y n )) with H(x) = 1 for x ≥ 0 and 0 else. A background rate r 0 is usually not considered within a linear CM. However, we introduced r 0 to avoid the immediate model falsification by outliers occurring in grid cells with zero or negative stress changes. In particular, we set r 0 = 1∕(N ⋅ K ⋅ t) , which means that we assume that one earthquake is associated with the background activity in the total time period and the reservoir volume. The Coulomb model is equivalent to the model of Dempsey and Suckale (2017) in the case that all faults are critically loaded. We call this model CM cr . To account for subcritical stress values in the beginning, we introduce the new parameter S 0 , which defines the initial stress deficit. Only after the stress increase exceeds S 0 in a given location, then faults are assumed to be critically stressed and the local seismicity rate starts to be proportional to the stress change, while it was zero before. This model is called CM subcr .
In all cases, we also test a modified version of the model which accounts for the fault density. In particular, we consider that the predicted rate is proportional to the normalized fault density factor, R w (t k , x n , y n ) = (x n , y n )R(t k , x n , y n ).
For comparison of the model fits, we use the Akaike Information Criterion, AIC = 2(M f − ln L) , to account for the different number of free parameters M f . The model with the smallest AIC value performs best.

Results
We calculated the results for different location errors = 0, 1, 3, 5, 7, 10 km and ∞ . In each case, we optimized the parameters of the seismicity models and compared the  Table 1. While the individual model parameters are discussed in "Discussion", we first compare the fits of the different models.
Best solutions are obtained, if location errors of 1 km for Poisson/CM cr and 3 km for CM subcr /RS are considered. Overall, the CM is found to describe the seismicity better than the Poisson model, but only performs well, if a stress deficit S 0 is considered ( CM subcr ). This is in agreement with the findings of Dempsey and Suckale (2017). If criticality is assumed from the beginning, the CM ( CM cr ) predicts the largest seismicity rates in the time span from 1975 to 1982, when no seismicity was observed but production was largest. The fit of the CM subcr is similarly good for both types of input data; performing slightly better for the pressure input. However, the RS model explains the observed seismicity best. The improvement of the RS model with respect to CM subcr is small for the compaction data, but significant for the pressure data. The best RS model based on pressure data yields a fit with a difference of AIC = 36 to the second best model. In all cases, besides model CM subcr based on compaction, the fits improve by considering fault density. They further improve when additionally the fault orientations are taken into account which yield the best results for the model classes. For the Poisson model for instance, considering the fault density improves the model by AIC = 32 . Taking the fault orientation into account improves the model further by AIC = 43 and considering a location error of 1 km gains AIC = 56 . Larger location error does not further improve the model.
In the following, we analyze the model fits in the time and space domain. Figure 5 shows the time series of the summed seismicity rates in the reservoir predicted by the models in comparison to the observed seismicity. Here the results for the models based on pressure changes (model A) are shown by solid lines, while the corresponding results for the compaction strain data are marked by dashed lines. For both sets of input data, the seismicity rates predicted by RS shows a delayed and smooth onset of seismicity with increasing rates over time, as well as the decrease of seismicity rate after reduction in production in 2014; capturing the main features of the observed seismicity.
Similarly, the CM subcr based on compaction data fits equally well. However, the seismicity rates predicted by CM subcr based on pressure changes predict a rather abrupt onset of activity and cannot fully reproduce the strong increase observed until 2014. Figure 6 shows spatial maps of the temporally integrated seismicity rates until 2014 for the CM subcr and the RS model based on pressure data. To illustrate the effect of the fault Table 1 Fit results of seismicity models for the complete field using Eq. (9)  In the case where the fault density is ignored (a, b), the layer thickness is visible in the seismicity map with a smooth gradient towards the northern area, where most events are predicted. Additionally, the effect of some major faults is already indirectly visible, because they act as barriers in the pressure field leading to rapid changes of predicted seismicity. But the influence of the layer thickness is weaker for the RS model in Fig. 6b, here the contrasts in pressure changes are stronger preserved, causing a slight damping of the predicted seismic density in the northern part of the field. For example at grid point ( X = 250 , Y = 610 ) 0.32 events per km 2 are predicted by CM subcr and RS predicts 0.27 events per km 2 . The lower two panels (c) and (d) of Fig. 6 show corresponding maps for the models with fault density (see "Seismicity model and testing approach" and Fig. 3a). The influence of the fault density is dominant for both models. The seismicity is focused to regions of high fault density resulting in a remarkable agreement to the locations of observed seismicity, which is found when the fault density weight is applied to the Poisson model. The seismicity map for the CM subcr model, however, overestimates the seismicity rates in the northwestern area ( X < 240 , Y > 600 ) with values over 0.25 events per km 2 , whereas the values of the RS model reach only 0.2 events per km 2 north of the reservoir outlines. In the RS model two main regions of seismicity are predicted, one in the centre of the reservoir and one in the SW (see Fig. 6d). In these regions most seismicity occurred. In detail, there are some misfits, e.g. for the SW region the seismicity is overestimated by the model. These could indicate heterogeneities in the field as suggested by Candela et al. (2019). Comparing the two input data sets (model A-pressure changes and B-compaction strain) the spatial pattern of pressure changes in the raw data do not match to the earthquake locations whereas the compaction data fit significantly better, although the maximum is offset relative to the region with highest seismicity (see Fig. 2). In general, the compaction data are more focused than the pressure changes. The compaction strain is high in the central area where the most seismicity occurs and in three areas in the southern part of the reservoir where only little seismicity is observed. These areas in the south are enhanced relative to the central area with the largest compaction because the reservoir layer is only between 50 m to 150 m thick compared to over 220 m resulting in higher strain values (see Figure S1). Additionally, here are high fault densities. The low seismicity may be related to thick salt layer above the reservoir in these regions. Furthermore, in the northern central area ( X < 245 , Y > 597 ) the predicted seismicity underestimates the observed seismicity despite there being some compaction and a high fault density. Because of these places of overestimated and underestimated seismicity, the compaction strain input yields worse results than the pressure data and the fault density weighting does not improve the results for compaction strain. Regarding the two models, the results of the RS model and the CM subcr model cannot be distinguished, as shown by the same AIC results (see Figure S1 and Table 1). Figure 7 shows four maps of successive time slices with predicted and observed seismicity to evaluate the spatial and temporal performance of the RS based on pressure data.
The four time slices represent distinct changes in the observed seismicity. The first time period is the longest from 1960-1991 and represents the period, when production was high but no seismicity was recorded (Fig. 7a). In this period, the predicted seismicity is very low ( < 0.02 events per km 2 ) although a period of 31 years is summed. The other three  Table 1. The M ≥ 1.5 earthquakes occurring in the corresponding time period are plotted as magnitude scaled turquoise circles. The corresponding results for the compaction data are shown in the Supplementary material (see Figure S3) time slices cover each approximately ten years for better comparison. In the course of the next 9 years from 1991 to 2000 (Fig. 7b), some widely distributed events with M ≤ 2.7 happened. The predicted seismicity shows zones of slightly elevated seismic density (up to 0.1 events per km 2 ) in several areas of the field matching with the locations of the recorded events. In the following time period from 2000 to 2010 (Fig. 7c), the modeled seismicity is more focused in some areas (up to 0.3 events per km 2 ) likewise the observed events are more spatially clustered. These areas are nearby but do not always perfectly match. In this time period five events with M ≥ 3.0 and one of these with M = 3.5 occurred which could have triggered some of the subsequent events with magnitudes M ≥ 1.5 , i.e. partly responsible for the observed offset. Elsewhere the modeled seismicity increases, around 0.1 events per km 2 in most parts of the reservoir. The time period from 2010 to 2017 includes the production change in 2014. The maximum predicted seismicity increases to a density of nearly 0.5 events per km 2 , expected in a slightly wider area extending more to the south. These areas are in general agreement with the locations of the observed events. The widening happens in the time period 2010 to 2014 for the predicted seismicity as well as for the observed seismicity, when more events occur in the southern region. After 2014 an overall decrease in predicted and observed seismicity is visible with slightly pronounced seismicity in the southwestern region of the field. The decease is partially caused by the shorter observation period of three years. In one region in the northern central area some earthquakes occurred, where the model predicts very low seismicity ( << 0.01 events per km 2 ). But if larger earthquake location uncertainties are considered, these events could fall in nearby areas, where large seismicity rates are predicted. Overall the predicted seismicity matches well to the temporal and spatial occurrence of the observed seismicity, despite the involved uncertainties (earthquake location error, model errors) and assumed homogeneity of seismicity parameters.
The significant difference between the AIC values of the CM subcr and RS model result from their different spatiotemporal predictions. The accordant four time slices for the CM subcr display these differences (see Supplement Figure  S2). Whereas the first time slice of the CM subcr captures the very low seismicity, the seismic density patterns for the time periods 1991-2000, 2000-2010 and 2010-2017 all resemble Fig. 7c and do not account for the observed changes of the seismicity pattern in the three time periods.
The results differ systematically for the RS model based on compaction strain (see Supplement Figure S3). Here little seismicity is expected for the time period until 1991, but some patches of elevated seismicity are already visible (up to 0.2 events per km 2 ), which grow in the following time periods to a larger size and a seismic density over 0.3 events per km 2 (2000: 0.3 events per km 2 ; 2010: 0.35 events per km 2 ; 2017: 0.32 events per km 2 ). For the time period 1991 to 2000, the regions do not correspond to the locations where the earthquakes occurred ( Figure S3b). In the following time slices ( Figure S3c and d), the model predicts high seismicity only for parts of the central region with most events and increasing seismicity in the southern region where only after 2010 some seismicity occurs ( Figure S3d). This is the reason why the model based on pore pressure changes has better fitting results. It appears the compaction strain model fails to capture the underlying mechanism of the induced seismicity.

Discussion
The modeling of the reservoir related seismicity at the Groningen gas field provides a unique opportunity to test different theoretical approaches and models and to identify which key parameters control the induced seismicity. However, it is important to check the model consistency and uncertainty. The parameter estimations might be inconsistent for different fitting periods and the parameter uncertainties can be large due to the rather small number of earthquake data. To address these points, we estimated the parameters for four fitting periods (1960-2000, 1960-2005, 1960-2010, 1960-2017), each time sampling all parameter sets leading to results comparable to the best fit in the same period. In particular, we chose those results with | AIC| ≤ 5 relative to the optimal parameters. Figure 8 shows the time series of the corresponding parameter sets for the different fitting periods in comparison to the observed evolution of the seismicity until 2019. To account for the aleatoric variability, the 90 % confidence interval is shown in gray for the case of the optimal parameter set. The confidence interval is calculated assuming a Poisson process with the mean value predicted by the RS model. The lower and upper boundary of the 90% confidence intervals of all curves with | AIC| ≤ 5 is additionally marked by the black lines. For the shortest fitting period until 2000 with 35 events (a) the confidence interval of the best result does not predict the observed seismicity rate, but the models with similar good fits span a wide range of rates, with larger rates that almost cover the observed trend until 2018. For the next fitting time period until 2005 with 67 events (b) the best results shift to larger seismicity rates and the parameter range is slightly larger than for the previous time period, but with later observed seismicity rates still underestimated. When the fitting time is extended to 2010 with 134 events the result changes (c). The uncertainties increase and the high seismicity rates of 2011 and 2013 are inside the spread of the curves with comparable fits until 2010. This evolution shows that in the time period between 2000 and 2010 the seismicity pattern changed significantly. One reason could be that some of the observed M ≥ 1.5 events are aftershocks, e.g. triggered by the mainshocks with M > 3 , which are not included in the model predictions. Bourne et al. (2018), e.g. found 10-20% aftershocks for the time span 2012-2017. Another reason could be that our assumption of uniform parameters in space is an oversimplification and the spatial migration of the seismic activity leads to other effective model parameter estimations.
For the fitting period until 2017, the parameters are rather well constrained and the discrimination between the alternative models is not easy. The observed seismicity is inside the 90 % confidence interval of the best result except for the four years 1994, 2006, 2011 and 2013. Such outliers are statistically expected because one out of ten data points should fall on average outside a 90% confidence interval. The results highlight the involved uncertainties and that for a limited data set, i.e. with 35 events until 2000, the full solution space needs to be considered and other solutions could be found. The parameter ranges are rather large for the full data set and are given by the following intervals, determined by the minimum and maximum values of all | AIC| ≤ 5-solutions of the spatiotemporal fit of the whole time period: A ∈ [ 1.0;2.0] in units of MPa/c 1 , t a ∈ [ 2.4 ⋅ 10 5 ; 2.0 ⋅ 10 9 ] years and r 0 ∈ [ 3 ⋅ 10 −8 ; 8 ⋅ 10 −4 ] events per year. The background rate is calculated for the whole volume of the reservoir, thus an average vertical extent of 225 m. Assuming a total seismogenic width of 10 km and a b value of 1.0, the tectonic background rate would be approximately 0.2 M ≥ 0 -events per year in the whole area, thus approximately only one M ≥ 0-event in 5 years. Therefore, the region would by typically be identified as aseismic, as often assumed. On the other hand, the effective normal stress at reservoir level can be roughly estimated by means of the difference between lithostatic and hydrostatic stress, which yields a value of approximately eff = 45 MPa at reservoir level. From lab experiments, the A value is known to be about 0.01 (Dieterich 1994). Thus, A = A eff ∕c 1 = 1.5 ± 0.5 leads to c 1 =Ã( S∕ h ) ∈ [ 0.2;0.5] according to Eq. (4). Assuming a Biot-Willis coefficient = 0.8 and a Poisson's ratio = 0.2 yields Ã = 0.6 (Lele et al. 2015) and thus S∕ h ∈ [ 0.4;0.8] , which is in a reasonable range. The tectonic stressing rate can be estimated as A eff ∕t a ≤ 1.9Pa∕yr . This low tectonic stressing rate in combination with the low background rate explains the very slow reaction of the fault system to stress changes. These results give us some confidence that the RS model captures the major physics and can be used to estimate future seismicity levels for hazard estimation. The improvement of the spatial fit by implementation of the fault density supports the model that the induced stresses accumulate at pre-existing faults. The additional improvement by implementation of the fault orientation relative to the tectonic stress regime supports that most of the events are normal faulting events in the reservoir layer.
In contrast, the modeled seismicity rates for the compaction data set clearly underestimate the observed seismicity for all fitting time periods before 2017 (see Figure (S3) with slowly increasing seismicity rates for the longer data sets. At the same time, the model parameter range is similarly narrow for all data subsets. Here the Fig. 8 Seismicity rates as a function of time. The predictions of the RS model based on pressure changes are shown for the best parameters and parameters with similar good results ( | AIC| ≤ 5 ). The time limit for the parameter fits is marked by the vertical line: a 1960-2000, b 1960-2005, c 1960-2010 and d 1960-2017. The 90% confidence interval related to the optimal fit is shown and additionally for the similarly good solutions displaying the uncertainties of the model results. The observed rate of M ≥ 1.5 earthquakes is shown as turquoise curve. The accordant results for the compaction data are shown in Supplement Figure S4 uncertainties of the model do not give a hint of the later evolution of the observations. This again implies that the underlying mechanism is not captured by the compaction data. Therefore, we do not discuss the resulting model parameters for the compaction model in detail. However, the values of t a and r 0 are comparable to model A and the differences in A are due to c 2 (see Eq. (6)) which depends on the dip angle of the faults and the seismic coupling coefficient. Information on fault dip and throw could improve the model with compaction strain but this goes beyond the scope of this work.
To give an impression of the future seismicity, Fig. 9 shows the prediction and the confidence interval of future seismicity based on the production scenarios of pressure changes provided by NAM. For comparison, the predicted seismicity rates for the compaction strain data are given in Figure S5 which show smaller rates than the predictions from a pressure change model.
The forecasts with pressure data predict an annual rate of M ≥ 1.5 earthquakes slowly decreasing to average values around 20 in 2024 but with a wide confidence range between 7 and 32. Depending on the production scenarios, the range shifts to smaller or higher rates. The realized production is close to the smallest scenario thus the low observed seismicity is reasonable.

Conclusion
No earthquakes were recorded in the Groningen gas field before gas production started. After a long delay relative to the start of production earthquakes started to occur, where earthquake frequencies roughly follow the temporal evolution of gas production rates. The M3.6 Huizinge earthquake in 2012 gave cause for serious concern and led eventually to the reduction of the gas production in 2014. For the estimation of the seismic hazard at the Groningen gas field, it is, therefore, crucial to have powerful datadriven models that allow to quantify the relation to the gas production.
In the past, the concept of rate-and-state (RS)-dependent friction has been used successfully in many case studies to model seismicity with a low number of parameters. For this reason, we designed a probabilistic model for the Groningen gas field, in which seismicity occurs as a statistical response on stress changes in a fault network with RS frictional behavior. As a main result, we find that all statistical features of observed seismicity in Groningen can be explained with our model that is governed by only three free parameters. This includes particularly the delayed and smooth occurrence of seismicity with increasing rates following the gas production rates. In this model, the delayed onset of observed seismicity is due to the very low background seismicity rate combined with the low tectonic stressing rate, corresponding to a long lasting frictional response to the fault system on stress disturbances. Using additional information on the fault density and orientation leads to an improvement of the fit. We also compare this model to the popular linear Coulomb failure model where the reservoir is initially subcritically loaded, i.e. the stress must first exceed a critical level. This model can explain the delayed occurrence of seismicity, but fails to describe the smooth increase of seismicity. Additionally, the model would impose that any tectonic stressing taking place would not be associated with a constant seismicity rate, whereas the RS results in a very low tectonic stressing rate which is in agreement with historical observations.
To summarize, the results suggest that the RS model is suitable to perform reasonable forecasts of future earthquake rates in Groningen in a simple and computationally inexpensive way.
We emphasize that our current model, only provides the frequencies of earthquakes. In future work, it would be interesting to extend the model towards the simulation of magnitudes. In a simple approach, random numbers from a given distribution, e.g. the Gutenberg-Richter distribution, can be drawn for each earthquake. This would allow for reevaluation of maximum expected earthquakes, as calculated in Zöller and Holschneider (2016) based on the RS model instead of the simple Poissonian process for earthquake occurrence.
Acknowledgements Open Access funding provided by Projekt DEAL. This work was funded by the Bundesministerium für Bildung und Forschung (BMBF) in the framework of the project SECURE Fig. 9 Predicted seismicity rates for the best model results based on pressure scenarios as a function of time for the three production scenarios between 2014 and 2024. The end of the fitting period for the parameter estimation is indicated by the vertical line. The model uncertainties are displayed like in Fig. 8. The accordant result for the compaction data is shown in Supplement Figure S5 (Sustainable dEployment and Conservation of Underground Reservoirs and Environment). We thank the Nederlandse Ardolie Maatschappij (NAM) for providing the data from the reservoir. The seismic catalog was kindly provided by Koninklijk Nederlands Meteorologisch Instituut (KNMI). GZ acknowledges support from the DFG Collaborative Research Center 1294 (project B04). We also gratefully acknowledge the suggestions of an anonymous reviewer that improved our initial manuscript and many thanks to Timothy Davis for the manuscript improvement. The figures were created with Generic Mapping Tools (GMT).
Funding The research was funded by BMBF (03G0872D) and supported by Deutsche Forschungsgemeinschaft (Grants SFB1294 project B04).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.