Original Article

Risk-based tolerable magnitudes for induced seismicity at subsurface operations in the Netherlands

Annemarie G. Muntendam-Bos^1,2 and Ryan Schultz³

¹Dutch State Supervision of Mines, The Hague, The Netherlands; ²Delft University of Technology, Delft, The Netherlands; ³Swiss Seismological Service, ETH Zürich, Zürich, Switzerland

Abstract

Induced seismicity is a risk that needs to be managed at any subsurface operation. Forward-thinking strategies should focus on the tolerable magnitude, M_TOL, which defines the level of nuisance, damage, and/or likelihood of fatality at which induced seismicity is no longer deemed acceptable. In this study, our objective is to derive risk-based M_TOL values for all active subsurface operations in the Netherlands, from producing gas fields to geothermal heat licenses and (CO₂) storage operations. In the Netherlands, the acceptability criteria generally adopted are for local risks: 1% probability of DS1 damage, 1% probability of DS2 damage, and a probability of an individual dying due to building collapse (Local Personal Risk) of 10^-5/year. We estimate M_TOL values based on these risk tolerances and find that M_TOL values are field/license specific. We compare our M_TOL estimates-based individual risk results to equivalent aggregate risks in terms of the number of households affected. This analysis shows that designing thresholds solely on local risk metrics would allow for a significant potential impact, in some cases even larger than the impact of the Huizinge earthquake, which arguably led to the abandonment of the large Groningen gas field. We propose that the combination of individual and aggregate risk acceptance criteria aimed at avoiding unacceptable damage and nuisance provides a socially more defendable approach to inform ‘red light’ thresholds in future traffic light protocols. These considerations will also be informative for guiding monitoring improvements, as they can highlight regions that are currently deficient.

Keywords: induced seismicity; tolerable magnitudes; gas production; geothermal; CO₂-storage

 

 

Introduction

Subsurface operations that perturb the local stress field carry the risk of inducing earthquakes. These include hydrocarbon production, geothermal energy, hydraulic fracturing, wastewater disposal, and natural gas storage. Comprehensive reviews of human-induced seismicity case histories have been published (Foulger et al., 2018; McGarr et al., 2002; Muntendam-Bos et al., 2022; National Research Council, 2013; Schultz et al., 2020b; Suckale, 2009). Ground motions sufficiently large to cause damage to nearby buildings and infrastructure occurred in several cases and have even led to regulatory intervention, for instance the significant limitations on injection rates imposed on wastewater disposal in Oklahoma (Grigoratos et al., 2021; Hincks et al., 2018), and in some cases, the termination of projects, such as the enhanced geothermal energy projects in Pohang, South Korea (Ellsworth et al., 2019), and Basel, Switzerland (Häring et al., 2008). However, even felt events alone have led to regulatory intervention. Examples include strict regulatory control on the hydraulic fracturing operations in Alberta, Canada (Alberta Energy Regulator (AER), 2015, 2019), the closing of the Castor gas storage project in Spain (Cesca et al., 2014, 2021), and the moratorium on shale gas extractions in the UK (Clarke et al., 2014; Kettlety et al., 2021).

In the Netherlands, induced seismic events have been observed at a variety of subsurface operations, including gas extraction, geothermal heat extraction, and underground gas storage (Muntendam-Bos et al., 2022). Though magnitudes have been very limited for underground gas storage and geothermal heat (magnitudes ≤ 1.7), felt and damaging events (up to magnitude 3.6) were induced during the extraction of hydrocarbons. A well-known example is the Groningen gas field, where multiple events between M_L = 3.0 and M_L = 3.6 resulted in extensive non-structural damage, a few cases of minor structural damage, significant social impact, and ultimately the cessation of gas production from the field (Muntendam-Bos et al., 2017, 2022; Van der Voort & Vanclay, 2015; Van Thienen-Visser & Breunese, 2015). But Groningen is not a stand-alone case; also the Bergermeer and Roswinkel gas fields have experienced events of comparable magnitudes and incurred significant damages (Muntendam-Bos et al., 2022).

Although, to date, magnitudes of events associated with gas storage and geothermal heat are small, the potential impact of being associated with damaging induced events on these emerging industries, which are vital for the transition to a more renewable energy system, emphasizes the importance of managing induced seismicity. Shutdown due to the occurrence of larger magnitude events, as occurred in other cases, would seriously hamper the growth of these two vital technologies.

Within the framework of seismic risk governance (Grigoratos et al., 2025; Zhou et al., 2024), preventive measures are taken to prevent/limit the extent of induced seismicity as much as possible, e.g. limiting the amount of stress induced on pre-existing, critically stressed faults through proper geological characterization of faults and limiting pressure and temperature changes due to injection/production. However, unknowns and uncertainties remain and require mitigation measures.

The de facto management of induced earthquakes is typically done through the implementation of a traffic light protocol (TLP) (Bommer et al., 2006; Schultz et al., 2021b; Verdon & Bommer, 2020). However, the effectiveness of TLPs depends on the tacit assumptions that induced earthquakes of unacceptable strength are preceded by small magnitude precursory events, and that the escalation to unacceptable magnitudes can be prevented by modifying or stopping subsurface operations (Baisch et al., 2019; Verdon & Bommer, 2020). Thus, TLP’s are inherently reactive; operational decisions (e.g. a reduction or cessation of injection/production rates) are taken in response to observed seismic events.

In a recent work, adaptive traffic light systems (ATLSs) have been proposed (Mignan et al., 2017). These systems use a forecasting model to estimate the probability of an unacceptable event occurring based on (near) real-time observed activity. Operations can be amended or stopped when this probability exceeds a given threshold. What these traffic light methods share is the need to properly derive pre-set magnitude or ground motion thresholds for ‘amber’ and ‘red’ lights based on some acceptable level of risk (Bommer et al., 2015; Grigoratos et al., 2025; Schultz et al., 2021a, 2021b, 2022b; Verdon et al., 2025).

Often, TLPs are magnitude-based due to its implementational simplicity and the fact that its derivation can be done in near real-time, e.g. Berlín, El Salvador (Bommer et al., 2006), Basel, Switzerland (Häring et al., 2008), Pohang, South Korea (Ellsworth et al., 2019), Helsinki, Finland (Ader et al., 2020), Bowland Shale, UK (Clarke et al., 2014), and hydraulic fracturing in British Columbia, Canada (Kao et al., 2018). In the Netherlands, the Dutch mining law stipulates that all subsurface activities located onshore and within 12 nautical miles (22 km and 224 m) from the Dutch coast (in this study further referred to as near-shore activities) are required to execute a seismic risk assessment and implement a seismic risk management system. In practice, a very basic, not risk-based, TLP is often part of this risk management system: The main onshore gas operator has implemented a TLP with pre-set uniform magnitude thresholds for all their gas fields and storages. For geothermal heat extraction, each operator translates policy-defined ground motion thresholds to site-specific magnitudes using only the expected value of the ground motion prediction equation (GMPE).

In this study, we adapt the risk-based approach of Schultz et al. (2022b) to better inform TLP-design for subsurface operations in the Netherlands. We define the magnitude at which the consequences become intolerable as M_TOL. For the TLP-design, the ‘red light’-magnitude (M_RED) is defined as the threshold requiring immediate suspension of operations. To account for the risk of subsequent trailing seismicity, M_RED should be set an appropriate level below M_TOL (Schultz et al., 2021b; Zhou et al., 2024). However, the former Ministry of Climate and Green Growth (MinCGG)¹ of the Netherlands has specified that they do not intend to require operators to account for trailing seismicity in the TLP-design. Therefore, we focus on estimating M_TOL and discussing the implications of this and other policy choices on the TLP-design.

Method

The risk-based TLP approach was first developed by Schultz et al. (2020a, 2021a, 2021b). In 2022, the method was adapted to build risk scenarios for possible future enhanced geothermal projects in the Netherlands (Schultz et al., 2022b). However, as we are not interested in incorporating trailing seismicity effects, we limit the implementation of the method to the approach applied in Verdon et al. (2025). Note that the determination of M_TOL is relatively insensitive to the operation type, so these results can be more broadly adapted. In addition, these M_TOL results could be used to inform operation-specific M_RED thresholds later.

The method (Figure 1) consists of three main components: (1) definition of the source and shake grids, (2) computation of the hazard posed by a hypothetical event at a given location, and (3) estimation of the possible impact of the given hazard. In a Monte Carlo approach, these components are combined, and risk scenarios are built for all near-shore and onshore producing gas fields, all Dutch gas production prospects, which are near-shore or onshore, all issued Dutch licenses for geothermal heat extraction, and all near-shore and onshore gas storage sites. Subsequently, we use the iso-risk approach to derive values of M_TOL (Schultz et al., 2020a, 2021a, 2021b, 2022b). Maps showing the locations of the gas fields, license areas, and storage sites are provided in Figures S1–S4 of the Supplementary Material.

Figure 1. Schematic overview of the workflow of our method to estimate M_TOL.

We note that our method does not incorporate any assumption on the probability of an earthquake of a particular magnitude occurring at a particular location. This assumption is justified for building a TLP, since we only need to know when to stop the operation. In the following subsections, each of the components is discussed in more detail. For further details, we refer the reader to previous work (Schultz et al., 2020a, 2021a, 2021b, 2022b; Verdon et al., 2025).

Source and shake grids

We sample possible earthquake source locations for each near- or onshore gas field, prospect, geothermal license area, and gas storage site. We start by defining a very dense ‘source grid’ of 0.0025° × 0.0025° (approximately 250 × 250 m) over the whole Netherlands and Dutch offshore. Subsequently, we retain only the grid locations within the contours of the specified fields or license areas. The very dense grid ensures that we include at least a single potential source location in all fields, even the very small gas fields. For computational purposes, the source grid of the larger fields and licenses is upscaled in a stepwise manner up to a maximum grid of 0.025° × 0.025° (approximately 2.5 × 2.5 km).

Besides computational arguments, the main reasoning for retaining only the grid points within the contours of the analyzed fields and licenses is the objective to include site-specific depth distributions for the earthquake sources. Detailed analysis of the high-precision hypocenter locations derived for Groningen events using full-waveform inversion (Willacy et al., 2019) and microseismic monitoring at reservoir depth in both the Groningen (Willacy et al., 2019) and the Bergermeer gas fields (Baisch et al., 2016) show that the induced events occur predominantly within the reservoir layer, at the top of the reservoir. Injection-induced events typically occur close to the injection interval (Buijze et al., 2019). Therefore, we adopt the top of the stratigraphic layer in which the gas (storage) reservoir or geothermal activity is located as the reference depth of the ‘source grid’ for that specific field/license.

It is not uncommon for gas or hot water to be produced from several, vertically stacked reservoirs. In general, the depth of the top of the stratigraphic unit of the shallowest reservoir layer is adopted as a reference depth. As shallower events are more impactful at the surface than deeper events, this is a conservative assumption. In some cases, production from a very thin (< 20 m) Vlieland sandstone layer (part of the late Jurasic, early Cretaceous Rijnland subgroup) is stacked with production from much thicker reservoir layers in the Permian Zechstein Carbonates and/or Rotliegend Sandstone formations. These specific shallow and thin reservoir layers are dominated by dip closures and have been found to be too geomechanically weak to induce felt seismic events (Muntendam-Bos, 2021; Van Eijs et al., 2006; Van Thienen-Visser et al., 2018). Therefore, for fields of this configuration, the top of the second reservoir layer is adopted to derive reservoir depth. For reference, a stratigraphic column for the Dutch subsurface is shown in Figure 2.

Figure 2. Stratigraphic column showing the geological sequence and groups specific for the Dutch subsurface. Note that for our main stratigraphic layers in the paper, we have combined the Nedersaksen and Schieland Late Jurassic Groups and use the abbreviation NS for this combination (Source: https://www.dinoloket.nl/stratigrafische-nomenclator).

By limiting the ‘source grid’ to the contours of the gas fields and storage sites, we assume that the stress changes induced by the activity, and which may subsequently result in seismic activity, are limited to these reservoirs. This means that stress changes induced by pressure and/or temperature changes in adjacent, hydraulically connected aquifers are not considered. Since little is known of the lateral extent of these aquifers or of the trailing diffusion behavior of the aquifer in response to the activity, it is our preference in this study to limit our analysis to the areas immediately, and most significantly affected by the subsurface operation: the reservoirs itself.

Apart from our exclusion of the very thin, shallow Vlieland gas bearing reservoir layers in the presence of much thicker deeper reservoir layers, our analysis does not explicitly incorporate the presence of seismogenic faults in our source areas. As noted earlier, the purpose of our analysis is to determine when an operation needs to stop, rather than if geological and operational conditions are encountered that could enable a M_TOL event.

Hazard calculation

In the second step, we calculate the ground motions of an earthquake with a given hypocenter and magnitude. We use the GMPE of Ruigrok et al. (2026). This empirical model, called JAPEK, is based on Dutch-induced events in the magnitude range of 1.0 ≤ ML ≤ 3.6, epicentral distances of 0–110 km, and nucleation depths of 2.3–4 km. The model contains empirical relations both for peak ground velocity (PGV), which is relevant for most of our risk metrics, and peak ground acceleration (PGA), which is relevant for the ‘likelihood of death’ risk metric.

The GMPE contains two correction terms. The first accounts for the significant additional damping of the ground motion amplitudes near the epicenter in the presence of the (thick) Zechstein salt layer between the hypocenter and the Earth’s surface. As a significant number of fields and licenses are located below or within the Zechstein group, we incorporate this correction based on stratigraphic unit for all reservoirs and licenses located in layers dating from the Carboniferous and Permian eras.

The second correction term accounts for the distance-dependent amplification in the top 200 m soil column. This term also accounts for the amplifying nature of the relatively weak shallow top-soils with an average Vs₃₀ of 200 m/s present in the Netherlands. The additional effect of local site amplification due to small deviations from this average Vs₃₀ (at the recorded locations ranging between 137 and 317 m/s) was investigated and found to have negligible effect on the performance of the models (Ruigrok et al., 2026). The explanation for this negligible effect was found in a clear magnitude dependence of the variability of the amplification, suggesting a strong frequency dependence of the local site response. This implies a much more complicated behavior of the site response than described by the simple Vs₃₀-proxy in the model. Given this much more complicated and frequency-dependent amplification behavior, which is not included in the model, in combination with the negligible effect on the model performance when implementing the Vs₃₀-proxy relation, made us decide to exclude an additional Vs₃₀-related local site amplification in our computation of the ground motions.

The likelihood of death or fatality chance (FC) is related to building collapse. Therefore, we need to consider a geometric average of the spectral acceleration periods, as it is sensitive to structural damage. In the absence of a GMPE for spectral acceleration, we convert PGA into spectral acceleration (Graizer & Kalkan, 2009). Subsequently, we perform a geometric averaging of the spectral acceleration periods over the periods 0.01s, 0.1s, 0.2s, 0.3s, 0.4s, 0.5s, 0.6s, 0.7s, 0.85s, and 1.0s, which range is consistent with the fatality risk studies for Groningen (Crowley & Pinho, 2020). As was shown in Schultz et al. (2022b), this conversion provides an adequate way to compute average spectral acceleration and is consistent with other work (Bommer & Alarcon, 2006; Crowley & Pinho, 2020).

Risk calculation

In the third step, the ground motions of step 2 are translated into risk estimates. We consider three risk metrics: nuisance, damage, and FC. We use the Community Decimal Intensity (CDI) to distinguish various levels of nuisance (Wald et al., 2012), where the levels 2–6 roughly correspond to the subjective descriptions ‘just felt’, ‘exciting’, ‘somewhat frightening’, ‘frightening’, and ‘extremely frightening’, respectively. Similarly, we consider different damage states (DSs) (Korswagen et al., 2019), where level 1 corresponds to visible light damage (0.1–1 mm cracks) and level 2 to easily observable light damage (1–2 mm cracks). Higher damage states (3–5) contain a level of structural damage and are only considered as part of the FC.

Nuisance/fragility functions relate the different levels of nuisance/damage to the level of ground motion. We use the nuisance functions developed by Schultz et al. (2021c), which were fitted to data from North America. There may be a systematic difference in the way ground motions are experienced between people in the central US and the people in the Netherlands. After all, the Netherlands is in a tectonically stable region with very little natural seismicity occurring only in the southeast of the country. This could introduce a systematic error in our nuisance risk estimation. However, we use these functions because these were developed on data from induced events in the Central and Eastern United States, which means the population is also not regularly exposed to, and therefore not used to, earthquake ground motions. Further work to develop functions specific to the experience of the Dutch population to earthquake-induced ground motions would be beneficial but is beyond the scope of this study.

For damage, we use the fragility functions of Korswagen et al. (2019), which are derived for induced seismicity in Groningen, the Netherlands. The Groningen building stock can be considered representative for the average building stock in the Netherlands. The fragility functions are representative of pre-1950 unreinforced masonry buildings. These buildings are among the most vulnerable types of building in the region. Thus, our approach ensures that the individual risk of these more vulnerable buildings determines M_TOL, and our results can be regarded as somewhat conservative for the average building stock.

The FC relates to the individual risk of a person dying because of buildings (partially) collapsing due to an earthquake. To compute FC, we use the predefined average-Groningen vulnerability function of Crowley and Pinho (2020), which relates FC to the average spectral acceleration.

To avoid implausible outcomes, we set the likelihood of nuisance to zero for PGV’s below 0.1 mm/s, the likelihood of damage to zero for PGV’s below 0.5 mm/s, and the likelihood of fatality to zero for average spectral accelerations below 0.05 g.

To compute aggregate risk metrics, we multiply the impact likelihood by the number of households at the given shake grid location. To estimate the number of households, we use the LandScan population model (Rose et al., 2019; Figure S5 of the Supplementary Material) with an average of three people per household. The total number of households affected is computed by summing over the affected number of households per grid point. For the individual risk metrics, we only use the population model to derive the nearest populated shake grid point and the epicentral distance to this point. Here, we note that for all near-shore fields and sites, we limit our analysis to populated grid cells onshore. Any possible damages to nearby infrastructure on the North Sea, such as wind farms, hydrocarbon pipelines, wells, or platforms, are currently not considered in this analysis and require separate risk assessments.

Monte Carlo sampling

Both the input parameters and the models of the components contain uncertainties and variabilities. We use a 1000-trial Monte Carlo sampling approach in which the input parameters and variabilities are randomly perturbed around their distributions (Figure 3). These distributions are described in more detail later.

Figure 3. Perturbations to the input parameters and model variabilities in our Monte Carlo analysis.

For each point on the source grid, the outcome from these repeated trials is a set of 1000 risk scenarios for each of our risk metrics, representing the statistical distributions of our metrics. We define M_TOL at a source grid location as the mean likelihood or mean number of impacted households at the given risk tolerance for that metric. By defining M_TOL based on the mean of the distributions instead of the median (50–50 chance) of a given risk, we focus on the commonly used and more informative expected amount of risk impact.

Earthquake depth variability

From the reference depth given by the top of the stratigraphic unit containing the gas (storage) reservoir or geothermal activity, the scenario depths are perturbed about a distribution that skews in favor of events just below the reference depth while also allowing for slightly deeper events. This distribution is justified by the observation from the Groningen gas field that events predominantly occur within and toward the top of the reservoir layer (Spetzler & Dost, 2017; Willacy et al., 2019), and injection-induced events typically occur close to the injection interval (Buijze et al., 2019).

Ground motion model uncertainty

The variability in the ground motion model accounts for both the within-event variability φ and the between-event variability τ. This yields a standard deviation , which is expressed in Ruigrok et al. (2026). The standard deviation for PGA is larger than for PGV. This stems from a larger between-event variability, which reflects a larger variability in the source term for higher frequencies. In the Monte Carlo iterations, the standard deviation is multiplied with the GMPE perturbation dGM.

Fragility function uncertainty

The nuisance model contains two terms dN₁ and dN₂, for which standard deviations are computed based on the fit to the observations used to calibrate the model. Perturbations are applied to both terms.

The fragility functions for damage include a pre-damage term Ψ0, which adapts the vulnerability of a building to new damage in the presence of already existing damage (Korswagen et al., 2019). We sample this term from a half-Gaussian distribution with a standard deviation of 0.15. This distribution represents pre-damage, which is imperceptible.

Finally, the uncertainty in the vulnerability function is considered via a Gaussian representation of the logic tree branches (Crowley & Pinho, 2020).

Acceptability of a given level of impact

The impact of an earthquake can be multifaceted and is generally categorized into types: such as the risk of a person dying, the damage to buildings and infrastructure, the level of nuisance experienced by the public, the financial impact from an induced event (Langenbruch et al., 2020), or a combination thereof. Except for the risk of a person dying as a result of an earthquake, which, in most risk policies, is considered as an individual risk metric (Marzocchi et al., 2015), acceptable levels are generally inferred in terms of aggregate numbers of households or buildings affected (Schultz et al., 2021a, 2022b; Verdon et al., 2025).

In the Netherlands, the Dutch risk policy and the Dutch mining law are more focused on individual risk criteria. For the Groningen gas field and geothermal heat extraction, the Dutch mining law defines a ‘safe’ tolerance for the annual risks of dying of a fictitious person due to an induced earthquake (known as Local Personal Risk, LPR) of 10^-5/yr. This risk tolerance is comparable to the levels defined by many organizations, which range between 10^-6 and 10^-4, depending on the source and consent to risk (Marzocchi et al., 2015). Though not formally stipulated in the mining law, this risk tolerance choice is generally also adopted for other subsurface operations such as gas extraction and storage. Therefore, we use the probability of 10^-5 as the tolerance level for our FC risk metric.

For damage, the Dutch mining law only stipulates that this should be avoided as much as reasonably possible. For geothermal heat extraction, the additional guidance is provided that some level of damage is allowed if the operator is insured for the financial impact of the damage incurred by an induced event. In the TLP-guideline for geothermal operators of MinCGG (https://www.nlog.nl/sites/default/files/2025-06/beleidskader_geothermie_v1306_0.pdf), this is translated to an ‘amber-light’ threshold at ground motions corresponding to a 1% probability of DS1 damage (3 mm/s; Roos et al., 2009; Staalduinen & Geurts, 1998; Staalduinen et al., 2019) and a ‘red light’ threshold at ground motions corresponding to a 1% probability of DS2 damage (33 mm/s; Korswagen et al., 2019). From this, we infer that for the damage-risk tolerance in the Netherlands, we probably should adopt a 1% probability of DS2 damage occurring due to an induced earthquake at any individual house. However, we argue that the intention of the ‘red light’ threshold of a TLP is to allow for regulatory intervention before damaging ground motions occur, and runaway events (events that grow much larger than the operator intended) only result in minor, non-structural damage. In the light of effective TLP-design, we also assess the more averse individual damage-risk tolerance of 1% probability of DS1 damage.

Even though the Dutch government focusses on individual risk metrics, our aim is to compute both individual and aggregate risk metrics. By computing both, we will be able to compare equivalent risks as well as gain insight into the impact of alternative approaches on M_TOL. Therefore, we also compute additional tolerances, where the number of households affected by either DS1 or DS2 damage exceeds 1. We will hereafter refer to this tolerance as [N(DS1) ≥ 1] and [N(DS2) ≥ 1], respectively.

Prior work on nuisance has led to empirically derived aggregate risk tolerances for Groningen of 33,556 households affected by CDI2, 19,844 households affected by CDI3, and 11,077 households affected by CDI4 (Schultz et al., 2022b). Compared to the tolerances derived similarly for the UK (Verdon et al., 2025) and for hydraulic fracturing in North America (Schultz et al., 2021b), these tolerances are on the high side, which may be due to the large population density in the Netherlands. Although some studies have shown that nuisance can be an important metric in the TLP-design (Cremen & Werner, 2020; Schultz et al., 2021a, 2022b), the Dutch government finds this approach too restrictive toward its objectives for the sustainable use of the Dutch subsurface in the energy transition. Therefore, we do not impose nuisance tolerances a priori. However, we will investigate the equivalent aggregate nuisance risk imposed by the tolerances of the other risk metrics in the discussion.

Results

We applied our approach to all near- and onshore producing gas fields, gas prospects, geothermal hear licenses, and gas storage sites. We first present the results for M_TOL obtained by imposing the individual risk acceptance criteria of 1% probability of DS1 damage, 1% probability of DS2 damage, and an FC of 10^-5. Subsequently, we will present the results for the aggregate risk acceptance criteria of [N(DS1) ≥ 1] and [N(DS2) ≥ 1].

Note that the results of fields/licenses/sites located shallower than 1 km depth are indicated in red in Tables S1–S8 of the supplementary Material as these are outside the usable depth range of the GMPE of 1–8 km (Ruigrok et al., 2026). At depths shallower than 1 km, surface waves are likely dominating the recordings, and surface-wave attenuation is not described by the GMPE model.

M_TOL based on individual risk acceptance criteria

In Figures 4–6, we plot the M_TOL values derived using the individual risk acceptance criteria for DS1 damage, DS2 damage, and FC for each of the gas fields, geothermal license areas, and gas storage sites. The values per field/license/site are also provided in Tables S1–S4 of the Supplementary Material. Here, we report the lowest M_TOL value derived on the source grid of each field/license/storage site as the objective is to derive a single M_TOL-value per risk metric for each field/license/storage site, which could be used in the TLP-design.

Figure 4. DS1 Iso-damage maps derived for a 1% probability of DS1 damage occurring at a building in the nearest populated shake grid point. See also Tables S1, S2, S3, and S4 in the Supplementary Material for the specific M_TOL values per field, license, and gas storage site.

 

Figure 5. DS2 Iso-damage maps derived for a 1% probability of DS2 damage occurring at a building in the nearest populated shake grid point. See also Tables S1, S2, S3, and S4 in the Supplementary Material for the specific M_TOL values per field, license, and gas storage site.

 

Figure 6. Iso-risk maps derived for a fatality probability of 10^-5 for a fictitious person residing in or near a building in the nearest populated shake grid point. See also Tables S1, S2, S3, and S4 in the Supplementary Material for the specific M_TOL values per field, license, and gas storage site.

The values derived for M_TOL based on the DS1 risk metric are generally about 0.7 magnitude units lower than those based on the DS2 risk metric. The difference between the M_TOL-values derived based on the FC risk metric and DS2 damage is relatively small, about 0.2 magnitude units. This small difference could be the result of the fact that the fragility functions relate the probability of damage to PGV, while the vulnerability functions use the PGA-derived average spectral accelerations. In the GMPE, the PGV and PGA show a different dependence on hypocentral distance and damping in the subsurface (Ruigrok et al., 2026). In addition, the fragility and vulnerability functions are non-linearly dependent on the ground motions. Especially at low probabilities, small changes in ground motion could significantly change the probability and hence the derived mean values of M_TOL.

The variability of M_TOL for a single risk metric tolerance shows a clear dependence on depth, the presence of the Zechstein salt above the reservoir, and distance to the shore (Figure 7). The additional damping of the salt layer results in larger values of M_TOL. The effect is most pronounced for the damage metrics, which are related to PGV. Again, this points to the different dependence of PGV and PGA on the attenuation of the subsurface layers finding its way into M_TOL. The influence of the distance to the nearest populated grid point becomes apparent in the scattering of the higher M_TOL values, clearly deviating from the depth-related trends. This is further confirmed by the dependence of M_TOL on hypocentral distance shown in Figure S6 of the Supplementary Material, which shows an alignment of these higher values with the observed trend.

Figure 7. The variability of M_TOL with depth for (a) DS1-damage, (b) DS2-damage, and (c) chance of fatality. The grey symbols indicate fields below the Zechstein salt layer, and the black symbols indicate fields located in layers above the Zechstein salt. The filled symbols correspond to fields with an epicentral distance of more than 2 km to the nearest populated grid point. PF: producing fields; PP: prospect fields; GT: geothermal heat licenses; SF: storage site locations.

M_TOL based on aggregate risk acceptance criteria

As for the individual risk acceptance criteria, the values of M_TOL per field/license/site are provided in Tables S5–S8 of the Supplementary Material. In Figure 8, we plot the M_TOL derived from the individual risk acceptance criteria versus those based on the aggregate risk acceptance criteria. In general, we find that M_TOL based on the aggregate risk acceptance criteria is lower than the values derived from the individual risk criteria by 0.1–0.4 magnitude units. Here, we recall that to compute aggregate risk metrics, we multiply the impact likelihood by the number of households at the given shake grid location and subsequently sum over the affected number of households per grid point. Hence, in the presence of a high number of households exposed, the aggregate [N(DS1) ≥ 1] / [N(DS2) ≥ 1] criteria are reached at lower magnitudes than the criteria of 1% chance of DS1/DS2 damage, respectively. For fields/sites offshore densely populated areas, such as the CO₂-storage sites offshore South-Holland near the city of The Hague, the low probabilities will aggregate even more quickly, and hence, the M_TOL values of the aggregate risk criteria are significantly lower than the corresponding magnitude based on the individual risk criteria, especially for DS1 damage (up to 0.7 magnitude units).

Figure 8. The M_TOL-values derived from the individual risk criteria versus those based on the aggregate risk acceptance criteria [N(DS1) ≥ 1] and [N(DS2) ≥ 1]. The black line indicates where the two estimates are the same. PF: producing fields; PP: prospect fields; GT: geothermal heat licenses; SF: storage site locations.

Only a couple of fields located on the transition from on- to offshore (e.g. the fields underneath the Wadden-islands in the north of the Netherlands) show lower values based on the individual risk criteria. This is most likely due to the relatively low population density on these islands, while the nearest populated grid cell is still immediately above the source location. The low number of households exposed leads to the criteria of 1% chance of DS1/DS2 damage being reached at lower magnitudes than the aggregate [N(DS1) ≥ 1] / [N(DS2) ≥ 1] criteria, respectively.

Discussion

We will now discuss the implications of our results for effective TLP-design.

Comparison between individual and aggregate risks

As described earlier, the Dutch risk policy and mining law are focused on individual risk acceptance criteria. Our results show that for most operations, the value of M_TOL based on the aggregate risk criteria is lower than the corresponding values based on the individual risk criteria. This immediately suggests that when the individual risk criteria for damage are adopted, some extent of damage to houses is to be expected. Furthermore, this suggests that aggregate damage metrics can provide complementary information for the TLP-design.

For a representative selection of gas fields and geothermal heat licenses, we have computed the number of households expected to experience damage when an earthquake of magnitude M_TOL(1% DS1) or M_TOL(1% DS2), as computed for this field or license, respectively, was to occur. The results are shown in Table 1. The number of households expected to experience damage clearly depends on the population density (see also Figure S5 of the Supplementary Material). An M_TOL event in the vicinity of major cities such as Groningen, Alkmaar, Delft, Leeuwarden, or The Hague leads to more potential households experiencing damage than an M_TOL event in a field or license located near smaller cities such as the gas fields Emmen, Roswinkel (both near the town of Emmen), and Eleveld (near the town of Assen), or in a more rural area such as the geothermal licenses Luttelgeest, Middenmeer, and Californië.

Table 1. Overview of the number of households (HH) expected to experience DS1 or DS2 damage when an earthquake corresponding to M_TOL(1% DS1) or M_TOL(1% DS2), respectively, was to occur.

Gas field	Depth (km)	SU	DS1 damage			DS2 damage
			MTOL (1% DS1)	#HH (1% DS1)	#HH (CDI3)	MTOL (1% DS2)	#HH (1% DS2)	#HH (1% DS1)	#HH (CDI3)
Bergermeer	–2.30	RO	2.7	22	23,618	3.4	15	4,262	89,839
Eleveld	–3.26	RO	2.9	13	15,666	3.6	12	2,888	72,217
Emmen	–3.14	ZE	2.8	3	11,160	3.5	15	2,618	44,099
Gaag	–3.02	RB	2.8	27	53,749	3.5	39	4,911	262,525
Groningen	–2.79	RO	2.8	8	28,216	3.5	28	6,727	77,658
Roswinkel	–2.11	RB	2.4	3	9,889	3.1	15	1,820	42,632
Geothermal license	Depth (km)	SU	DS1 damage			DS2 damage
			MTOL (1% DS1)	#HH (1% DS1)	#HH (CDI3)	MTOL (1% DS2)	#HH (1% DS2)	#HH (1% DS1)	#HH (CDI3)
Californië-IV	–1.78	DI	2.3	2	3,656	3.0	2	257	27,113
Delft I	–1.94	KN	2.3	134	36,924	3.1	162	10,635	146,552
Klazienaveen	–2.86	RB	2.6	15	6,588	3.3	29	1,419	29,998
Kwintsheul	–2.07	NS	2.4	27	32,169	3.1	47	4,274	143,222
Leeuwarden	–2.45	RO	2.8	9	16,167	3.5	22	3,909	50,637
Luttelgeest	–1.80	RO	2.7	4	3,388	3.3	4	417	24,045
Middenmeer-I	–2.17	RO	2.7	2	7,348	3.4	2	752	42,823
The representative gas fields and licenses have been selected based on depth variation as well as variability in population density in the vicinity of the operation. The abbreviations of the stratigraphic units (SU) correspond to the abbreviations of the group names in Figure 2. #HH(1% DS1) – number of households expected to experience DS1 damage; #HH(DS2) – number of households expected to experience DS2 damage; #HH(CDI3) – number of households expected to feel the earthquake (CDI3 nuisance).

We recall that the current regulations for geothermal TLP-design define M_RED at the 1% chance of DS2 damage. Therefore, for our representative selection of fields and licenses, we have also computed the number of households expected to experience DS1 damage when an earthquake of M_TOL (1% DS2) was to occur (#HH(1% DS1) in Table 1). In most cases, the number of households affected by DS1 damage falls in the range of 1,000–10,000 households. These values are an order of magnitude larger than the number of households impacted by DS1 damage computed at ‘red light’ thresholds for various hydraulic fracturing sites in North America (Schultz et al., 2021a, 2021b) and well above the provisional damage tolerance of T_DS1 = 183 derived by Schultz et al. (2022b) for the Netherlands based on the analysis of the Huizinge and Roswinkel events. Given a representative average damage compensation for the repair of slight, non-structural DS1-damage of 6,000–9,000 Euros per household (derived from experience in the Groningen earthquake area; Instituut Mijnbouwschade Groningen, 2026), geothermal operators are currently liable for potential damage claims ranging from 6M€ to 90M€. It is highly questionable whether a geothermal heat operator has the financial resources to be liable for these prospective damage claims and/or insurance to cover these costs can be obtained.

The equivalent nuisance impact of this M_TOL choice (#HH(CDI3) in Table 1) is even more extensive. The number of households that are expected to feel the M_TOL(1% DS2) event ranges from tens of thousands to more than 250,000. Schultz et al. (2022b) also derived provisional nuisance tolerances for the Netherlands, which, for CDI3, amounted to 19,844 households. The nuisance impacts of all our selected fields and licenses are well above this tolerance. When a 1% chance of DS1 damage would be imposed for M_RED, the equivalent number of households derived is spread around the nuisance tolerance value. In urban areas, the equivalent nuisance impact is still exceeding the tolerance, while in rural areas, the impact is well below the tolerance.

Comparison with historical earthquake impact

Table 2 provides an overview of all historical, induced M ≥ 3.0 earthquakes in the Netherlands. For some of these events, the number of recognized damage claims is known and included in the overview. For eight of these, we compute the consequences of an event of this magnitude at its specific hypocenter location. For this calculation, we used the observed ground motions to calibrate the ground motion variability dGM (Figure S7 of the Supplementary Material) We note that the number of ground motion data points used for the calibration was very limited for the events up to the 2012 Huizinge earthquake (ranging from < 5 for the Roswinkel and Westeremden events to 8 for Huizinge). For the later events, the number of data points exceed 50, which is an immediate result of the significant extension of the Groningen monitoring network in the period 2014–2016. All other steps in our procedure were kept identical, and again, we used 1000 Monte Carlo iterations to derive the posterior uncertainty distribution. The results are shown in Figure 9.

Table 2. Catalog of all observed induced earthquakes M_L ≥ 3.0 in the Netherlands.

Date	Location	ML	PGVmax (mm/s)	Latitude	Longitude	Depth (km)	# Recognized damage claims
2025-11-14	Zeerijp	3.4	35.1	53.348	6.774	3	2,117‡*
2022-10-08	Wirdum	3.1	9.2	53.319	6.772	3
2021-11-16	Garrelsweer	3.2	15.5	53.309	6.751	3
2019-05-22	Westerwijtwerd	3.4	8.3	53.328	6.652	3	10,001†
2018-01-08	Zeerijp	3.4	28.1	53.363	6.751	3	7,389†
2015-09-30	Hellum	3.1	9.8	53.234	6.834	3	6,100†
2014-02-13	‘t Zandt	3	14.4	53.3568	6.7817	3
2013-07-02	Garrelsweer	3	5.8	53.2935	6.785	3
2013-02-07	Zandeweer	3.2	14.5	53.3892	6.6667	3
2012-08-16	Huizinge	3.6	34.5	53.3453	6.6717	3	1,800*
2011-06-27	Garrelsweer	3.2	12.4	53.3028	6.7867	3
2009-05-08	Zeerijp	3	4.9	53.3538	6.7617	3
2008-10-30	Loppersum	3.2	14.4	53.3367	6.72	3
2006-08-08	Westeremden	3.5	12.5	53.3503	6.6967	3	275
2003-11-10	Stedum	3		53.3253	6.69	3	43
2003-10-24	Garrelsweer	3		53.295	6.7917	3	5
2001-09-10	Bergen	3.2		52.6526	4.7118	2.5	316*
2001-09-09	Bergen	3.5		52.651	4.7133	2.5
2000-10-25	Roswinkel	3.2	31	52.8318	7.0517	2.2
1998-07-14	Roswinkel	3.3	46	52.8325	7.0533	2.2	141
1997-02-19	Roswinkel	3.4	60	52.8323	7.0384	2.2	204
1994-09-21	Bergen	3.2		52.658	4.708	2.5	120
1994-08-06	Bergen	3		52.6543	4.711	2.5	54
The number of recognized damage claims has been derived from de Crook et al. (1998), Dost & Kraaijpoel (2013), KNMI (1994a, 1994b), Postmes et al. (2020), and Roos et al. (2009), and peak ground velocities from ground motion databases (Dost & Haak, 2007; Dost & Kraaijpoel, 2013; Ntinalexis et al., 2022). †Since the Huizinge earthquake, there is an increased awareness of earthquakes and damage caused by induced earthquakes. *For these events, a few cases of DS2 damage have been reported and recognized. ‡Number of damage reports received by the Institute of Mining Damage Groningen (IMG) on 18-11-2025.

 

Figure 9. Earthquake impact scenarios. Comparison of DS1 damage impact of (a) selected events M ≥ 3.0 before the Huizinge 2012-08-16 event and (b) selected M ≥ 3.0 events after the 2012-08-16 Huizinge event. Vertical dashed lines indicate the corresponding (color of dashed line and bar diagram) number of recognized damage claims. Events have been selected based on (1) availability of ground motion information to calibrate the ground motion variability and (2) availability of damage claim information.

Schultz et al. (2022b) identified the 2012-08-16 Huizinge event as the most impactful and the archetype of an intolerable event. It was also the event that triggered social perception of seismic risk and policy changes (Van der Voort & Vanclay, 2015). Therefore, we contrast the events prior to this event (Figure 9a) against the M ≥ 3.0 events after the Huizinge event (Figure 9b), where the Huizinge event serves as a reference. Of the older earthquakes (Figure 9a), the number of recognized damage claims of the Roswinkel 1997-02-19 and Roswinkel 1998-07-14 events is significantly lower than predicted by the model (well below the 10th percentile). The predicted distributions coincide with the Huizinge distribution, but significantly less damage reports were submitted and recognized (Table 2). The reason for this could be that (induced) earthquakes and earthquake damage were still a relatively unfamiliar phenomenon in the Netherlands in the 1990s. For the Westeremden 2006-08-08 and Huizinge 2012-08-16 events, the number of recognized damage claims is higher than the expected values predicted by the model (41 and 676, respectively) and around the 90th percentile of the distributions, corresponding to 305 and 1717, respectively. The significant difference in damage impact between the Westeremden, which, up to the Huizinge event, was the strongest event observed in the Groningen gas field, and the Huizinge event corroborates the public perception of the Huizinge event being notably different from everything they had experienced before.

The number of recognized damage claims of the post-Huizinge earthquakes (Figure 9b) is much higher than predicted by the model and well outside the predicted distributions. Postmes et al. (2020) studied the strong increase in damage claims of the events following the Huizinge earthquake and concluded that earthquake awareness is awakened by the first damaging event occurring in an area. As awareness is heightened, damage of subsequent events gets reported quicker, and the extent of the reported damage is generally more extensive. In addition, they observed that with the increased seismic activity in the Groningen gas field, the number of claims attributable to a single event gets more diffuse. For instance, the number of damage claims of the 2019-05-22 Westerwijtwerd event was most likely blended with damage claims of the 2019-06-09 M = 2.5 earthquake at Garrelsweer. We note that the computed impact of the recent 2025-11-14 Zeerijp event is comparable to the impact of the Huizinge event. This is due to the quite significant ground motions associated with this event, especially in comparison to the 2006-08-08 Westeremden and 2019-05-22 Westerwijtwerd events.

The earthquake impact scenario comparison demonstrates that our model generally underestimates the number of DS1 damage reports, especially in regions where awareness of the possible occurrence of induced seismicity and associated damage has been raised. This further indicates that adopting the 1% chance of DS2 damage individual risk metric acceptance criteria as the ‘red light’ risk metric in the TLP-design leaves the operator vulnerable to very significant potential financial liabilities, which could be well beyond the operators’ financial resources.

Perspectives on TLP-design

So far, we have focused our analysis on M_TOL and the equivalent aggregate implications. The values of M_TOL represent the tipping point at which the seismic activity becomes unacceptable. In this section, we will discuss how M_TOL can inform effective TLP-design and address different perspectives on which risk acceptance criteria to adopt. Within this section, we discuss five inter-related topics: we reiterate the framework for ‘good practice’ red light design, provide critiques of the current Dutch TLP regulations, critiques of using individual risk for red light definitions, and address the implications for monitoring.

Framework for ‘good practice’ ‘red light’ design

Within the framework of a TLP, operations can proceed as normal during a ‘green light’. Once the threshold for the ‘amber light’ has been exceeded, mitigation strategies should be implemented with the aim to avoid further escalation of the seismic activity. After a ‘red light’, operations must be suspended indefinitely. In this simple context, adopting M_TOL as M_RED would seem reasonable. However, induced earthquakes experience jumps in magnitude of the next largest event and may persist after the suspension of operations. For the Dutch gas fields, trailing seismicity is regularly observed and directly related to the equilibration of reservoir pressures. Trailing event magnitudes are generally of the same order as the magnitudes observed during the operation period (see Table S10 of the Supplementary Material). However, in many instances of injection-induced seismicity, the largest and most impactful events occurred after suspension of the operation (Ellsworth et al., 2019; Langenbruch et al., 2011; Langenbruch & Shapiro, 2010; Schultz et al., 2017, 2023; Verdon & Bommer, 2020). It is therefore recommendable to set the value of M_RED at a value lower than M_TOL, for instance, by including a representative statistical trailing seismicity model (Schultz et al., 2022a). For operations that lack adequate observational data, this may entail physical modeling constraints.

Dutch TLP regulations

Following Dutch law and policy, this paper has mainly focused on risk from the perspective of individual risk exposure. The most risk-based TLP currently in practice in the Netherlands is the generic TLP for geothermal heat operations, which places M_RED at a ground motion loosely connected to a 1% chance of DS2 damage and M_AMBER at a ground motion loosely connected to a 1% chance of DS1 damage (Ministerie van Klimaat en Groene Groei, 2026). The underlying assumption is that by placing M_RED at this ground motion, the probability of ground motions which could possibly pose an actual fatality risk is sufficiently low. As this assumption is not informed by any quantitative risk assessment, we verify its validity against our risk-based results. Our results show that for all cases assessed, M_TOL(FC) exceeds M_TOL(1% DS2). This would suggest that indeed, a TLP using M_TOL(1% DS2) as M_RED could be sufficient to abide to the fatality norm. However, as indicated in the introduction, trailing seismicity has not been included in our analysis. If the increase in magnitude of the likely trailing events were to exceed +0.2 magnitude units, M_TOL(FC) would be exceeded, and a probability of incompliance exists. This probability of incompliance will predominantly be controlled by the amount of trailing seismicity an operation encounters and will range from insignificant to problematic. As a demonstrative example, in cases with a b-value of 1.0 and 90-10 split between operational/trailing events, the largest magnitude event would be +0.2, +1.0, and +2.0 units larger than M_RED with a chance of 67, 11, and 1%, respectively. For reference, the average ratio of operational/trailing events for the no longer producing Dutch near-/onshore gas fields is a 92-8 split, whereas the ratio for the Groningen gas field is a 95-5 split.

We recall that regarding the acceptability of damage, the Dutch law only states that it must be avoided as much as possible, which, for geothermal heat operations, is appended by the condition that the operator needs to have sufficient financial resources or insurance to compensate the predicted possible damages for which it is liable. Our equivalent aggregate impact results show that in the case of M_TOL(1% DS2) for M_RED, these liabilities could be quite extensive, especially in urban areas, and it is highly questionable whether operators, indeed, can bear this liability. Furthermore, the significant equivalent aggregate nuisance impacts of these moderate earthquakes to the regional community proximal to the subsurface operations could seriously jeopardize the social license to operate.

We argue that the intention of the ‘red light’ TLP-threshold is to provide a window for regulatory intervention before unacceptable risk is incurred. With this intention, ‘red light’ thresholds should be set such that the occurrence of events resulting in damage is rare and limited to non-structural, cosmetic damage only. Following the Dutch regulators’ preference for individual risk metrics, setting a ‘red light’ threshold at 1% chance of DS1 damage (M_TOL(1% DS1)) would then be more appropriate than the use of M_TOL(1% DS2). This alternative threshold does not eliminate the occurrence of damage, especially since trailing seismicity is not considered. On the other hand, the number of households predicted to incur damage is limited and of non-structural, cosmetic nature only, so the operator should be able to bear the financial liability or obtain an insurance to cover it.

At the same time, the equivalent nuisance impact is spread around the derived tolerance for the Netherlands (Schultz et al., 2022b). As expected, the nuisance impact is higher in urban areas and much lower in rural areas. Here, we additionally note that the nuisance functions were derived for events of M ≥ 3.1 in Central US (parts of Kansas, Oklahoma, and Texas; Schultz et al., 2021c), and the authors concluded that due to the impulsive nature of the smaller induced events (M < 3.9), more vigorous ground shaking is required to increase the person’s confidence the event is felt. In urban areas, the low levels of ground motion associated with the M_TOL (1% DS1) magnitudes could therefore be lost in the increased background noise levels. If we offset this to the fact that the nuisance tolerance was derived for the rural Groningen region, this may not be representative for the nuisance tolerance of the densely populated Dutch urban areas. Hence, there may be some justification for setting a ‘red light’ threshold at 1% chance of DS1 damage (M_TOL (1% DS1)) and accepting the calculated higher nuisance impact in the urban areas. Alternatively, it may be worthwhile to explore individual nuisance at various tolerance levels, for example M_TOL (50% CDI3) or M_TOL (10% CDI4). This way nuisance is not completely ignored while also not being overly restrictive.

Similarly, the ‘amber light’ (M_AMBER) defines the low end of the magnitude range in which operators intervene and apply mitigation strategies to avoid escalation to unacceptable levels. Hence, M_AMBER should be set at an adequate magnitude difference below the ‘red light’ threshold M_RED. Statistical considerations of induced seismicity sequences (Schultz et al., 2020a) and observations of large numbers of sequences in practice (Verdon & Bommer, 2020; Watkins et al., 2023) lead to a general recommendation of a magnitude difference of approximately two magnitude units between M_AMBER and M_RED. This band is necessary to minimize the probability of a run-away event to jump straight from ‘green’ to ‘red’. As stated earlier, the current Dutch TLP regulations place M_AMBER at a ground motion loosely connected to a 1% chance of DS1 damage (Ministerie van Klimaat en Groene Groei, 2026). Our results show that the differences between the M_TOL(1% DS1) values and the M_TOL(1% DS2) values are about 0.7 magnitude units. This difference is insufficient to allow for proper, timely mitigation once exceeding M_AMBER to avoid escalation toward magnitudes exceeding M_RED.

The use of individual or aggregate risk metrics for TLP-design

Alternatively, adopting an aggregate risk metric acceptance criterion would further differentiate TLP-design between operators by including the local population density as a factor to be considered. The main advantage of this approach would be that the financial liability of all operators becomes identical, and the social license to operate of the different industries is best secured. The equivalent aggregate nuisance impact for the selected fields associated with the [N(DS1) ≥ 1] acceptance criteria (Table S9 of the Supplementary Material) shows that only for the Gaag gas field, located in the very densely populated Rotterdam-The Hague area, the number of households impacted just exceeds the nuisance tolerance (#HH(CDI3) = 20,359 versus T_CDI3 = 19,844). The drawback of this approach is that the acceptance criteria is not sufficiently low to also ensure that in low-populated areas, the risk of the individual persons remains acceptable. Our results show that in these low-populated areas, M_TOL based on the individual risk metric is lower than M_TOL based on the aggregate metric. An iso-Risk combination map taking the minimum ‘red light’ value of the respective individual and aggregate risk acceptance criteria for each field would accommodate the advantages of both risk metrics. Figure S8 of the Supplementary Material shows the iso-Risk combination map indicating for each field, license, and storage site, which DS1 acceptance criteria (individual or aggregate) result in the lowest M_TOL value.

Implications for monitoring

Assuming that the DS1 iso-Risk-criteria combination map for M_RED was adopted, then values of M_RED would be as low as M = 1.6. In combination with the appropriate difference of two magnitude units, this implies that M_AMBER should be set as low as M = –0.4. The Dutch jurisdictions may opt to loosen the constraint to some extent, but for most fields/licenses, monitoring down to a completeness of at least M = 0.6 would be required for a properly functioning TLP. In the Netherlands, the current magnitude of completeness of the network is at least M = 2.0 country wide. In most regions with subsurface activities, the completeness of the network is at least M = 1.5, regionally going down to M = 1.0 and more locally to M = 0.5 (see Figure S9 of the Supplementary Material). The Dutch Seismological Survey (KNMI) is currently in the planning stage of extending the network such that the network of completeness everywhere is at least M = 1.5 (except for the far southwest and some Wadden-islands). However, even after the extension operators would still need to rigorously assess the completeness of the network at their specific field or license and, if necessary, impose local improvements essential to properly manage the seismic risk.

Further considerations

The approach described in this paper is entirely from a regulatory perspective. No operational information has informed the TLP-design, nor does it consider any real-time seismic information. As operations progress, near real-time cataloguing of induced events provides a wealth of information. Analysis of spatiotemporal variations in event rates and Gutenberg-Richter b-values can help delineate fault structures, provide insight into the seismogenic nature of the reservoir, and inform predictions of the next largest credible earthquake. These factors can be used in adaptive seismic risk management and inform mitigation strategies (Schultz et al., 2024; Zhou et al., 2024).

Furthermore, when selecting criteria for decision-making, it is important to clarify whether the priority is given to individual safety or overall societal safety. Future studies could benefit from incorporating public consensus and stakeholder engagement, involving authorities, industry operators, and the public in jointly defining acceptable risk criteria and priorities.

Conclusions

In this study, we computed maximum tolerable magnitudes (M_TOL) for all Dutch near- and onshore gas fields, prospects, geothermal heat licenses, and gas storage locations. We adapted the risk-based method and implementation of Schultz et al. (2022b) that was developed for the Netherlands. Particularly, we replaced the Groningen-specific GMPE of Bommer et al. (2022) by the newly developed GMPE of Ruigrok et al. (2026), which is specific to the Dutch subsurface. The derived M_TOL magnitudes can inform a TLP, which not only is tailored to the Dutch risk policy but also ensures the social license to operate is protected. Our results show that under the current risk policy, the aggregate nuisance and damage impacts are sufficiently large to potentially jeopardize the social license to operate of future operations. We propose that the combination of individual and aggregate risk acceptance criteria aimed at avoiding unacceptable damage and nuisance provides a socially more defendable approach to inform ‘red light’ thresholds in future TLPs.

Author declaration

AMB – Conceptualization, Formal analysis, Software, Funding acquisition, Investigation, Visualization, Writing – original draft; RS – Methodology, Software, Writing – review & editing.

The authors declare the following financial interests/personal relationships that may be considered as potential competing interests: Annemarie Muntendam-Bos is employed by the Dutch State Supervision of Mines (SodM), the Dutch safety authority for subsurface operations, and acts as expert on issues pertaining induced seismicity. SodM has not had any input on the conception, development, or analysis presented in this study. The expressed opinions on the effective TLP-design are those of the authors and do explicitly not represent the formal position of SodM. Ryan Schultz has no conflicts of interest to declare.

Acknowledgments

We would like to thank an anonymous reviewer and the editor for providing valuable comments that helped improve this manuscript. All models and input information were derived from earlier studies and public sources: GMPE model (Ruigrok et al., 2026), population density, nuisance functions (Schultz et al., 2021c), fragility function (Korswagen et al., 2019), vulnerability function (Crowley & Pinho, 2020), earthquake catalog (www.knmi.nl/kennis-en-datacentrum/dataset/aardbevingscatalogus), ground motion databases (Dost & Haak, 2007; Dost & Kraaijpoel, 2013; Ntinalexis et al., 2022), and information on gas fields, geothermal licenses, and gas storage sites as well as stratigraphic maps (www.nlog.nl). Routines and publicly available data used to create the figures and results of this paper are available online at GitHub (https://github.com/amuntendambos/LTEnl). Additional materials on the workflow are also discussed in prior studies (Schultz et al., 2021a, 2021b).

Supplementary materials

Supplementary material associated with this article can be found in the online version.

References

Ader, T., Chendorain, M., Free, M., Saarno, T., Heikkinen, P., Malin, P.E., Leary, P., Kwiatek, G., Dresen, G., Bluemle, F. & Vuorinen, T., 2020. Design and implementation of a traffic light system for deep geothermal well stimulation in Finland. Journal of Seismology 24: 991–1014. https://doi.org/10.1007/s10950-019-09853-y

Alberta Energy Regulator (AER), 2015. Subsurface order no. 2: monitoring and reporting of seismicity in the vicinity of hydraulic fracturing operations in the Duvernay zone. Fox Creek.

Alberta Energy Regulator (AER), 2019. Subsurface order no. 7: monitoring and reporting of seismicity in the vicinity of hydraulic fracturing operations in the Duvernay zone. Red Deer.

Baisch, S., Koch, C. & Muntendam-Bos, A., 2019. Traffic light systems: to what extent can induced seismicity be controlled? Seismological Research Letters 90: 1145–1154. https://doi.org/10.1785/0220180337

Baisch, S., Koch, C., Vörös, R. & Rothert, E., 2016. Induced seismicity in the bergermeer field: hypocenter relocation and interpretation. Report TAQA003. Bad Bergzabern: QCon.

Bommer, J.J. & Alarcon, J.E., 2006. The prediction and use of peak ground velocity. Journal of Earthquake Engineering 10: 1–31. https://doi.org/10.1080/13632460609350586

Bommer, J.J., Crowley, H. & Pinho, R., 2015. A risk-mitigation approach to the management of induced seismicity. Journal of Seismology 19: 623–646. https://doi.org/10.1007/s10950-015-9478-z

Bommer, J.J., Oates, S., Cepeda, J.M., Lindholm, C., Bird, J., Torres, R., Marroquín, G. & Rivas, J., 2006. Control of hazard due to seismicity induced by a hot fractured rock geothermal project. Engineering Geology 83: 287–306. https://doi.org/10.1016/j.enggeo.2005.11.002

Bommer, J.J., Stafford, P.J., Ruigrok, E., Rodriguez-Marek, A., Ntinalexis, M., Kruiver, P.P., Edwards, B., Dost, B. & van Elk, J., 2022. Ground-motion prediction models for induced earthquakes in the Groningen gas field, the Netherlands. Journal of Seismology 26: 1157–1184. https://doi.org/10.1007/s10950-022-10120-w

Buijze, L., Van Bijsterveldt, L., Cremer, H., Paap, B., Veldkamp, H., Wassing, B.B.T., Van Wees, J.-D., Van Yperen, G.C.N., Ter Heege, J.H. & Jaarsma, B., 2019. Review of induced seismicity in geothermal systems worldwide and implications for geothermal systems in the Netherlands. Netherlands Journal of Geosciences 98: e13. https://doi.org/10.1017/njg.2019.6

Cesca, S., Grigoli, F., Heimann, S., González, Á., Buforn, E., Maghsoudi, S., Blanch, E. & Dahm, T., 2014. The 2013 September–October seismic sequence offshore Spain: a case of seismicity triggered by gas injection? Geophysical Journal International 198: 941–953. https://doi.org/10.1093/gji/ggu172

Cesca, S., Stich, D., Grigoli, F., Vuan, A., López-Comino, J.Á., Niemz, P., Blanch, E., Dahm, T. & Ellsworth, W.L., 2021. Seismicity at the Castor gas reservoir driven by pore pressure diffusion and asperities loading. Nature Communications 12: 4783. https://doi.org/10.1038/s41467-021-24949-1

Clarke, H., Eisner, L., Styles, P. & Turner, P., 2014. Felt seismicity associated with shale gas hydraulic fracturing: the first documented example in Europe. Geophysical Research Letters 41: 8308–8314. https://doi.org/10.1002/2014GL062047

Cremen, G. & Werner, M.J., 2020. A novel approach to assessing nuisance risk from seismicity induced by UK shale gas development, with implications for future policy design. Natural Hazards and Earth System Sciences 20: 2701–2719. https://doi.org/10.5194/nhess-2020-95

Crowley, H. & Pinho, R., 2020. Report on the fragility and consequence models for the Groningen field (Version 7). NAM-report, Nederlandse Aardolie Maatschappij (Assen).

de Crook, T., Haak, H.W. & Dost, B., 1998. Seismisch Risico in Noord-Nederland. KNMI-report TR-205. De Bilt: Royal Netherlands Meteorological Institute (KNMI).

Dost, B. & Haak, H.W., 2007. Natural and induced seismicity. In: T.E. Wong, D.A.J. Batjes & J. de Jager (eds.). Geology of the Netherlands. Amsterdam: Royal Netherlands Academy of Arts and Sciences, pp. 223–239.

Dost, B. & Kraaijpoel, D., 2013. The August 16, 2012 earthquake near Huizinge (Groningen). KNMI-report. De Bilt: Royal Netherlands Meteorological Institute (KNMI).

Ellsworth, W.L., Giardini, D., Townend, J., Ge, S. & Shimamoto, T., 2019. Triggering of the Pohang, Korea, earthquake (Mw 5.5) by enhanced geothermal system stimulation. Seismological Research Letters 90(5): 1844–1858. https://doi.org/10.1785/0220190102

Foulger, G.R., Wilson, M.P., Gluyas, J.G., Julian, B.R. & Davies, R.J., 2018. Global review of human-induced earthquakes. Earth-Science Reviews 178: 438–514. https://doi.org/10.1016/j.earscirev.2017.07.008

Graizer, V. & Kalkan, E., 2009. Prediction of spectral acceleration response ordinates based on PGA attenuation. Earthquake Spectra 25: 39–69. https://doi.org/10.1193/1.3043904

Grigoratos, I., Bazzurro, P., Rathje, E. & Savvaidis, A., 2021. Time-dependent seismic hazard and risk due to wastewater injection in Oklahoma. Earthquake Spectra 37: 2084–2106. https://doi.org/10.1177/8755293020988020

Grigoratos, I., Schultz, R., Van Ginkel, J., Gunatilake, T., Wiemer, S., Van Der Wal, J.L.N. & Muntendam-Bos, A.G., 2025. A generic seismic risk protocol for energy production sites. Bulletin of Earthquake Engineering 23: 1325–1347. https://doi.org/10.1007/s10518-024-02088-4

Häring, M.O., Schanz, U., Ladner, F. & Dyer, B.C., 2008. Characterisation of the Basel 1 enhanced geothermal system. Geothermics 37: 469–495. https://doi.org/10.1016/j.geothermics.2008.06.002

Hincks, T., Aspinall, W., Cooke, R. & Gernon, T., 2018. Oklahoma’s induced seismicity strongly linked to wastewater injection depth. Science 359: 1251–1255. https://doi.org/10.1126/science.aap7911

Instituut Mijnbouwschade Groningen (IMG), 2026. Schadevergoeding voor Groningers van 2012 tot nu. Accessed 21 January 2026. https://www.schadedoormijnbouw.nl/nieuws/2022/03/schadevergoeding-voor-groningers-van-2012-tot-nu

Kao, H., Visser, R., Smith, B. & Venables, S., 2018. Performance assessment of the induced seismicity traffic light protocol for northeastern British Columbia and western Alberta. The Leading Edge 37: 117–126. https://doi.org/10.1190/tle37020117.1

Kettlety, T., Verdon, J.P., Butcher, A., Hampson, M. & Craddock, L., 2021. High-resolution imaging of the ML 2.9 August 2019 earthquake in Lancashire, United Kingdom, induced by hydraulic fracturing during Preston New Road PNR-2 operations. Seismological Research Letters 92: 151–169. https://doi.org/10.1785/0220200187

KNMI, 1994a. Seismische analyse van de aardbeving bij Alkmaar op 21 september 1994. KNMI-report TR-167. De Bilt: Royal Netherlands Meteorological Institute (KNMI).

KNMI, 1994b. Seismische analyse van de aardbeving bij Alkmaar op 6 augustus 1994. KNMI-report TR-166. De Bilt: Royal Netherlands Meteorological Institute (KNMI).

Korswagen, P.A., Longo, M., Meulman, E. & Rots, J.G., 2019. Crack initiation and propagation in unreinforced masonry specimens subjected to repeated in-plane loading during light damage. Bulletin of Earthquake Engineering 17: 4651–4687. https://doi.org/10.1007/s10518-018-00553-5

Langenbruch, C., Dinske, C. & Shapiro, S.A., 2011. Inter event times of fluid induced earthquakes suggest their Poisson nature: inter event times of fluid induced earthquakes. Geophysical Research Letters 38: L21302. https://doi.org/10.1029/2011GL049474

Langenbruch, C., Ellsworth, W.L., Woo, J. & Wald, D.J., 2020. Value at induced risk: injection-induced seismic risk from low-probability, high-impact events. Geophysical Research Letters 47: e2019GL085878. https://doi.org/10.1029/2019GL085878

Langenbruch, C. & Shapiro, S.A., 2010. Decay rate of fluid-induced seismicity after termination of reservoir stimulations. Geophysics 75: MA53–MA62. https://doi.org/10.1190/1.3506005

Marzocchi, W., Iervolino, I., Giorgio, M. & Falcone, G., 2015. When is the probability of a large earthquake too small? Seismological Research Letters 86: 1674–1678. https://doi.org/10.1785/0220150129

McGarr, A., Simpson, D. & Seeber, L., 2002. 40 case histories of induced and triggered seismicity. W.H.K. Lee, H. Kanamori, P.C. Jennings, C. Kisslinger (eds.). International Handbook of Earthquake and Engineering Seismology, Part A. Elsevier, pp. 647–661. https://doi.org/10.1016/S0074-6142(02)80243-1

Mignan, A., Broccardo, M., Wiemer, S. & Giardini, D., 2017. Induced seismicity closed-form traffic light system for actuarial decision-making during deep fluid injections. Scientific Reports 7: 13607. https://doi.org/10.1038/s41598-017-13585-9

Ministerie van Klimaat en Groene Groei (minKGG), 2026. Beleidskader geothermie. Accessed 21 January 2026. https://www.nlog.nl/sites/default/files/2025-06/beleidskader_geothermie_v1306_0.pdf

Muntendam-Bos, A.G., 2021. Geomechanical characteristics of gas depletion induced seismicity in the Netherlands. In: Proceedings of the 55th US Rock Mechanics/Geomechanics Symposium, ARMA 21-1038. Presented at the 55th US Rock Mechanics/Geomechanics Symposium, Houston, TX, 20–23 June 2021.

Muntendam-Bos, A.G., Hoedeman, G., Polychronopoulou, K., Draganov, D., Weemstra, C., Van Der Zee, W., Bakker, R.R. & Roest, H., 2022. An overview of induced seismicity in the Netherlands. Netherlands Journal of Geosciences 101: e1. https://doi.org/10.1017/njg.2021.14

Muntendam-Bos, A.G., Roest, J.P.A. & De Waal, H.A., 2017. The effect of imposed production measures on gas extraction induced seismic risk. Netherlands Journal of Geosciences 96: s271–s278. https://doi.org/10.1017/njg.2017.29

National Research Council, 2013. Induced seismicity potential in energy technologies. Washington, DC: National Academies Press. https://doi.org/10.17226/13355

Ntinalexis, M., Kruiver, P.P., Bommer, J.J., Ruigrok, E., Rodriguez-Marek, A., Edwards, B., Pinho, R., Spetzler, J., Hernandez, E.O., Pefkos, M., Bahrampouri, M., van Onselen, E.P., Dost, B. & van Elk, J., 2022. A database of ground motion recordings, site profiles, and amplification factors from the Groningen gas field in the Netherlands. Earthquake Spectra 39(1): 687–701. https://doi.org/10.1177/87552930221140926

Postmes, T., de Jong, M., Schoutens, A.M.M., Adams, W.J., van Delden, C.N., Boendermaker, M.M., Kanis, B. & Stroebe, K., 2020. Door schade wijzer: Het proces van schademeldingen en haar geschiedenis. Groningen: RUG en Gronings Perspectief.

Roos, W., Waarts, P. & Wassing, B.B.T., 2009. Kalibratiestudie schade door aardbevingen. TNO report TNO-034-DTM-2009-04435.

Rose, A.N., McKee, J., Urban, M.L., Bright, E.A. & Sims, K.M., 2019. LandScan 2018 high-resolution global population data set. Oak Ridge, TN: Oak Ridge National Laboratory (ORNL).

Ruigrok, E., Kruiver, P.P., Van Thienen-Visser, K., Muntendam-Bos, A.G., Van Der Wal, J.L.N. & Vossen, C.E.J., 2026. Ground motion prediction equations for upper-crustal seismicity and soft sediments. Journal of Seismology 30: 31. https://doi.org/10.1007/s10950-026-10384-6

Schultz, R., Beroza, G., Ellsworth, W. & Baker, J., 2020a. Risk-informed recommendations for managing hydraulic fracturing–induced seismicity via traffic light protocols. Bulletin of the Seismological Society of America 110: 2411–2422. https://doi.org/10.1785/0120200016

Schultz, R., Beroza, G.C. & Ellsworth, W.L., 2021a. A risk-based approach for managing hydraulic fracturing–induced seismicity. Science 372: 504–507. https://doi.org/10.1126/science.abg5451

Schultz, R., Beroza, G.C. & Ellsworth, W.L., 2021b. A strategy for choosing red-light thresholds to manage hydraulic fracturing induced seismicity in North America. JGR Solid Earth 126(12): e2021JB022340. https://doi.org/10.1029/2021JB022340

Schultz, R., Ellsworth, W.L. & Beroza, G.C., 2022a. Statistical bounds on how induced seismicity stops. Scientific Reports 12: 1184. https://doi.org/10.1038/s41598-022-05216-9

Schultz, R., Ellsworth, W.L. & Beroza, G.C., 2023. An ensemble approach to characterizing trailing-induced seismicity. Seismological Research Letters 94: 699–707. https://doi.org/10.1785/0220220352

Schultz, R., Muntendam-Bos, A., Zhou, W., Beroza, G.C. & Ellsworth, W.L., 2022b. Induced seismicity red-light thresholds for enhanced geothermal prospects in the Netherlands. Geothermics 106: 102580. https://doi.org/10.1016/j.geothermics.2022.102580

Schultz, R., Quitoriano, V., Wald, D.J. & Beroza, G.C., 2021c. Quantifying nuisance ground motion thresholds for induced earthquakes. Earthquake Spectra 37(2): 789–802. https://doi.org/10.1177/8755293020988025

Schultz, R., Rinaldi, A.P., Roth, P., Madritsch, H., Gunatilake, T. & Wiemer, S., 2024. Pre-screening of induced seismicity risks for CO2 injection at Trüllikon, Switzerland. International Journal of Greenhouse Gas Control 138: 104239. https://doi.org/10.1016/j.ijggc.2024.104239

Schultz, R., Skoumal, R.J., Brudzinski, M.R., Eaton, D., Baptie, B. & Ellsworth, W., 2020b. Hydraulic fracturing-induced seismicity. Reviews of Geophysics 58: e2019RG000695. https://doi.org/10.1029/2019RG000695

Schultz, R., Wang, R., Gu, Y.J., Haug, K. & Atkinson, G., 2017. A seismological overview of the induced earthquakes in the Duvernay play near Fox Creek, Alberta. JGR Solid Earth 122: 492–505. https://doi.org/10.1002/2016JB013570

Spetzler, J. & Dost, B., 2017. Hypocenter estimation of induced earthquakes in Groningen. Geophysical Journal International 209(1): 453–465. https://doi.org/10.1093/gji/ggx020

Staalduinen, P. & Geurts, C., 1998. De relatie tussen schade aan gebouwen en lichte, ondiepe aardbevingen in Nederland: inventarisatie. TNO rapport 97-CON-R1523-1.

Staalduinen, P., Rots, J. & Terwel, K., 2019. Onderzoek naar de oorzaken van bouwkundige schade in Groningen Methodologie en case studies ter duiding van de oorzaken. TU Delft Report – Revisie 2nd of August of 2019. Delft University of Technology (Delft).

Suckale, J., 2009. Induced seismicity in hydrocarbon fields. In: R. Dmowska (ed.). Advances in geophysics. Elsevier, pp. 55–106. https://doi.org/10.1016/S0065-2687(09)05107-3

Van der Voort, N. & Vanclay, F., 2015. Social impacts of earthquakes caused by gas extraction in the Province of Groningen, the Netherlands. Environmental Impact Assessment Review 50: 1–15. https://doi.org/10.1016/j.eiar.2014.08.008

Van Eijs, R.M.H.E., Mulders, F.M.M., Nepveu, M., Kenter, C.J. & Scheffers, B.C., 2006. Correlation between hydrocarbon reservoir properties and induced seismicity in the Netherlands. Engineering Geology 84: 99–111. https://doi.org/10.1016/j.enggeo.2006.01.002

Van Thienen-Visser, K. & Breunese, J.N., 2015. Induced seismicity of the Groningen gas field: history and recent developments. The Leading Edge 34(6): 664–671. https://doi.org/10.1190/tle34060664.1

Van Thienen-Visser, K., Roholl, J.A., van Kempen, B.M.M. & Muntendam-Bos, A.G., 2018. Categorizing seismic risk for the onshore gas fields in the Netherlands. Engineering Geology 237: 198–207. https://doi.org/10.1016/j.enggeo.2018.02.004

Verdon, J.P. & Bommer, J.J., 2020. Green, yellow, red, or out of the blue? An assessment of Traffic Light Schemes to mitigate the impact of hydraulic fracturing-induced seismicity. Journal of Seismology 25: 301–326. https://doi.org/10.1007/s10950-020-09966-9

Verdon, J.P., Schultz, R. & Edwards, B., 2025. Tolerable magnitudes for induced seismicity at offshore carbon capture and storage projects. International Journal of Greenhouse Gas Control 142: 104335. https://doi.org/10.1016/j.ijggc.2025.104335

Wald, D.J., Quitoriano, V., Worden, C.B., Hopper, M. & Dewey, J.W., 2012. USGS ‘did you feel it?’ Internet-based macroseismic intensity maps. Annals of Geophysics 54: 5. https://doi.org/10.4401/ag-5354

Watkins, T.J.M., Verdon, J.P. & Rodríguez-Pradilla, G., 2023. The temporal evolution of induced seismicity sequences generated by low-pressure, long-term fluid injection. Journal of Seismology 27: 243–259. https://doi.org/10.1007/s10950-023-10141-z

Willacy, C., Van Dedem, E., Minisini, S., Li, J., Blokland, J.-W., Das, I. & Droujinine, A., 2019. Full-waveform event location and moment tensor inversion for induced seismicity. Geophysics 84: KS39–KS57. https://doi.org/10.1190/geo2018-0212.1

Zhou, W., Lanza, F., Grigoratos, I., Schultz, R., Cousse, J., Trutnevyte, E., Muntendam-Bos, A. & Wiemer, S., 2024. Managing induced seismicity risks from enhanced geothermal systems: a good practice guideline. Reviews of Geophysics 62: e2024RG000849. https://doi.org/10.1029/2024RG000849

Footnote

1 In 2026, the Ministry of Climate and Green Growth was merged into the Ministry of Economic Affairs and Climate.