We present an overview of induced seismicity due to subsurface engineering in the Netherlands. Our overview includes events induced by gas extraction, underground gas storage, geothermal heat extraction, salt solution mining and post-mining water ingress. Compared to natural seismicity, induced events are usually small (magnitudes ≤ 4.0). However, due to the soft topsoils in combination with shallow hypocentres, in the Netherlands events exceeding magnitude 1.5–2.0 may be felt by the public. These events can potentially damage houses and infrastructure, and undermine public acceptance. Felt events were induced by gas production in the north of the Netherlands and by post-mining water ingress in the south-east. Notorious examples are the earthquakes induced by gas production from the large Groningen gas field with magnitudes up to 3.6. Here, extensive non-structural damage incurred and public support was revoked. As a consequence, production will be terminated in 2022 leaving approximately 800 billion cubic metres of gas unexploited. The magnitudes of the events observed at underground gas storage, geothermal heat production and salt solution mining projects have so far been very limited (magnitudes ≤ 1.7). However, in the future larger events cannot be excluded. Project- or industry-specific risk governance protocols, extensive gathering of subsurface data and adequate seismic monitoring are therefore essential to allow sustainable use of the Dutch subsurface now and over the decades to come.

The Netherlands is in the midst of an energy transition with hydrocarbon production gradually declining, whereas the role of sustainable energy technologies is on the rise. One of these technologies is geothermal energy production from porous reservoirs at 1.5–3 km depth. As the number of geothermal projects increases, there is a growing concern that felt and/or damaging induced seismic events could occur as a result of geothermal operations. Over the last two decades, such events have occurred in the Netherlands due to gas production, notably in the Groningen gas field. However, the occurrence of felt events is limited to hydrocarbon fields in certain regions or reservoirs. Understanding where and for which plays these events are observed helps to estimate seismogenic potential for geothermal operations and other sustainable subsurface activities. Here, we summarise and review the main similarities and differences in terms of geological and geomechanical characteristics between the hydrocarbon and geothermal plays in the Netherlands, and we consider the differences in pressure and temperature changes. By doing so, we provide better insights into the factors that could play a role for fault reactivation and induced seismicity, and how these differ for hydrocarbon production and geothermal operations in the Netherlands. The review shows that geological characteristics for most geothermal target reservoirs are similar to those of hydrocarbon, albeit geothermal projects so far target higher porosity rocks than hydrocarbon reservoirs. On the other hand, pressure and temperature changes are very different, with significant depletion for hydrocarbon fields vs significant cooling around geothermal injection wells. The different operations result not only in different expected stress change magnitudes but also in a distinct spatio-temporal stress build-up on faults, which has implications for seismogenic potential and monitoring of these different operations.

We reanalysed the induced seismicity data from the Groningen gas reservoir. We used the well-maintained induced event catalogue of the KNMI. The distributions of seismic moments and interevent times show a power law behaviour over several decades, and we find that upon increasing the magnitude threshold, these distributions remained scale-invariant. Because of this scale-invariance, we can put a constraint on the average loading of the elastic energy within the reservoir, which upon reaching a critical value gives rise to the seismic events. We find that the elastic energy roughly increases proportional to time. We also propose a new machine learning approach for declustering the seismic events, separating correlated from independent events. We find that only few events are truly independent, i.e. exponentially distributed over time. There are also few aftershocks following an Omori-type power law. The bulk of events presents a Gamma distribution for interevent times. This gives us an indication that the rigidity in the reservoir is high, but whether this results in overall correlated events should be settled with physics-based arguments rather than statistical ones.

The occurrence of felt earthquakes due to gas production in Groningen has initiated numerous studies and model attempts to understand and quantify induced seismicity in this region. The whole bandwidth of available models spans the range from fully deterministic models to purely empirical and stochastic models. In this article, we summarise the most important model approaches, describing their main achievements and limitations. In addition, we discuss remaining open questions and potential future directions of development.

We address aseismic fault slip and the onset of seismicity resulting from depletion-induced or injection-induced stresses in reservoirs with pre-existing vertical or inclined faults. Building on classic results, we discuss semi-analytical modelling techniques for fault slip including dislocation theory, Cauchy-type singular integral equations and the use of Chebyshev polynomials for their solution and an eigenvalue-based stability analysis. We consider slip patch development during depletion for faults with zero, constant static and slip-weakening friction, and our results confirm earlier findings based on numerical simulation, in particular the aseismic growth of two slip patches that may subsequently merge and/or become unstable resulting in nucleation of seismic slip. New findings include improved approximate expressions for the induced seismic moment per unit strike length and a description of the effect of coupling between the slip patches which affects both forward simulation and eigenvalue computation for high values of the ratio of fault throw to reservoir height. Our implementation based on analytical inversion and semi-analytical integration with Chebyshev polynomials is more efficient and more robust than our numerical integration approach. It is not yet well suited for Monte Carlo simulation, which typically requires sub-second simulation times, but with some further development that option seems to be within reach. Moreover, our results offer a possibility for embedded fault modelling in large-scale numerical simulation tools.

Two closely spaced geothermal doublets were operated in the Californië geothermal field near Venlo, the Netherlands. The geothermal wells target the Dinantian Zeeland formation below 2 km depth. For several years, hot fluid was produced from the Tegelen fault, a regional fault in the Roer Valley rift system, until a felt M1.7 earthquake led to the suspension of geothermal activities. The Californië showcase provides a rare opportunity to retrospectively evaluate the assessment and the management of induced seismicity risks for a geothermal project. A seismic hazard assessment was conducted at several stages of the project, and seismicity was continuously monitored with a local station network.

In this paper, we report on the characteristics of the induced seismicity and evaluate the findings of the seismic hazard assessments conducted prior to the earthquakes. Seismic hazard assessments were based on numerical simulations of subsurface stress changes associated with geothermal operations. A geomechanical analysis indicated that the mapped faults in the subsurface are likely to be critically stressed. The largest hazard was inferred to result from thermo-elastic stresses, originating from cold water injection close to the Tegelen fault.

Subsequent earthquakes predominantly occurred near a production well after stopping or reducing production. We attributed this observation to a thermo-elastic stress load caused by cold water injection close to the Tegelen fault, combined with a counter-acting stabilisation of the fault due to pressure depletion during production. This mechanism was consistent with the dominating mechanism considered in the preceeding seismic hazard assessments. Although geothermal operations have not resumed yet, the geomechanical analysis indicates that re-locating one of the injection wells further away from the Tegelen fault could provide an efficient measure for mitigating induced seismicity risks at Californië.

Long-term exploration of the Groningen gas field in the Netherlands led to induced seismicity. Over the past nine years, an increasingly sophisticated Ground Motion Model (GMM) has been developed to assess the site response and the related seismic hazard. The GMM output strongly depends on the shear-wave velocity (VS), among other input parameters. To date, VS model data from soil profiles (Kruiver et al., Bulletin of Earthquake Engineering, 15(9): 3555–3580, 2017; Netherlands Journal of Geosciences, 96(5): s215–s233, 2017) have been used in the GMM. Recently, new VS profiles above the Groningen gas field were constructed using ambient noise surface wave tomography. These so-called field VS data, even though spatially limited, provide an independent source of VS to check whether the level of spatial variability in the GMM is sufficient. Here, we compared amplification factors (AF) for two sites (Borgsweer and Loppersum) calculated with the model VS and the field VS (Chmiel et al., Geophysical Journal International, 218(3), 1781–1795, 2019 and new data). Our AF results over periods relevant for seismic risk (0.01–1.0 s) show that model and field VS profiles agree within the uncertainty range generally accepted in geo-engineering. In addition, we compared modelled spectral accelerations using either field VS or model VS in Loppersum to the recordings of an earthquake that occurred during the monitoring period (ML 3.4 Zeerijp on 8 January 2018). The modelled spectral accelerations at the surface for both field VS and model VS are coherent with the earthquake data for the resonance periods representative of most buildings in Groningen (T = 0.2 and 0.3 s). These results confirm that the currently used VS model in the GMM captures spatial variability in the site response and represents reliable input for the site response calculations.

We perform statistical analyses on spatiotemporal patterns in the magnitude distribution of induced earthquakes in the Groningen natural gas field. The seismic catalogue contains 336 earthquakes with (local) magnitudes above $1.45$ , observed in the period between 1 January 1995 and 1 January 2022. An exploratory moving-window analysis of maximum-likelihood b-values in both time and space does not reveal any significant variation in time, but does reveal a spatial variation that exceeds the $0.05$ significance level.

In search for improved understanding of the observed spatial variations in physical terms we test five physical reservoir properties as possible b-value predictors. The predictors include two static (spatial, time-independent) properties: the reservoir layer thickness, and the topographic gradient (a measure of the degree of faulting intensity in the reservoir); and three dynamic (spatiotemporal, time-dependent) properties: the pressure drop due to gas extraction, the resulting reservoir compaction, and a measure for the resulting induced stress. The latter property is the one that is currently used in the seismic source models that feed into the state-of-the-art hazard and risk assessment.

We assess the predictive capabilities of the five properties by statistical evaluation of both moving window analysis, and maximum-likelihood parameter estimation for a number of simple functional forms that express the b-value as a function of the predictor. We find significant linear trends of the b-value for both topographic gradient and induced stress, but even more pronouncedly for reservoir thickness. Also for the moving window analysis and the step function fit, the reservoir thickness provides the most significant results.

We conclude that reservoir thickness is a strong predictor for spatial b-value variations in the Groningen field. We propose to develop a forecasting model for Groningen magnitude distributions conditioned on reservoir thickness, to be used alongside, or as a replacement, for the current models conditioned on induced stress.