Induced seismicity risk analysis of the hydraulic stimulation of a geothermal well on Geldinganes, Iceland

The rapid increase in energy demand in the city of Reykjavik has posed the need for an additional supply of deep geothermal energy. The deep-hydraulic (re-)stimulation of well RV-43 on the peninsula of Geldinganes (north of Reykjavik) is an essential component of the plan implemented by Reykjavik Energy to meet this energy target. Hydraulic stimulation is often associated with fluid-induced seismicity, most of which is not felt on the surface but which, in rare cases, can be a nuisance to the population and even damage the nearby building stock. This study presents a first-of-its-kind pre-drilling probabilistic induced seismic hazard and risk analysis for the site of interest. Specifically, we provide probabilistic estimates of peak ground acceleration, European microseismicity intensity, probability of light damage (damage risk), and individual risk. The results of the risk assessment indicate that the individual risk within a radius of 2 km around the injection point is below 0.1 micromorts, and damage risk is below 10⁻², for the total duration of the project. However, these results are affected by several orders of magnitude of variability due to the deep uncertainties present at all levels of the analysis, indicating a critical need in updating this risk assessment with in situ data collected during the stimulation. Therefore, it is important to stress that this a priori study represents a baseline model and starting point to be updated and refined after the start of the project.

1 Introduction

The city of Reykjavik, the capital and center of population of Iceland, meets 99.9 % of its district heating demand by geothermal energy (Gunnlaugsson et al., 2000). However, the growing population and the booming number of tourists are pushing the current supply of energy to its limit, since no new low-temperature wells have been drilled since 2001. In particular, additional sources of low-temperature heat need to be accessed to ensure a reliable heat provision for the city center. Therefore, there is an urgent need to increase the current capacity by drilling new low-temperature wells and stimulating older inactive wells.

One potential area for new low-temperature geothermal field developments is Geldinganes. Geldinganes is a peninsula within the city limits of Reykjavik (Fig. 1). The exceptional geothermal gradient in this area triggered the drilling of a well (RV-43) in 2001 after a gabbro body was identified as potential heat source and drilling target for this deviated well. Despite the required temperatures being reached, the flow rates were insufficient for economic production. At present, Reykjavik Energy (Orkuveita Reykjavíkur – OR) has reassessed this field for development of geothermal energy with new production wells. To additionally enhance the production, it is foreseen to hydraulically re-stimulate well RV-43 in order to improve its productivity to economical levels. In particular, a three-staged cyclic pulse stimulation is planned that will last for (circa) 12 d. The stimulation is expected to enhance productivity in three pre-existing fracture zones penetrated by RV-43 and isolated with straddle packers. Packer technology, commonly used in the oil and gas industry, is expected to be employed for the upcoming stimulation. Straddle packers not only allow isolating and injecting in selected narrow zones but also allow adjusting the straddle distance between the upper and the lower injection points.

Figure 1 Map view of the Geldinganes peninsula, the injection well, and seismic network. Source of the map: map data © 2019 Google Maps.

Like all energy technologies, the exploitation of deep geothermal energy is not risk-free. Therefore, an essential part of the implementation and licensing is a quantitative risk assessment comparable to existing regulations for health, safety, and environment (HSE) procedures. This analysis allows balancing the (perceived and real) risks against the (perceived and real) benefits. Over the last decade, induced seismicity has emerged as one of the risks – and often the most dominant one – to be faced (Giardini, 2009; Grigoli et al., 2017) in implementing industrial underground technologies (e.g., geothermal energy exploitation, water impoundment, CO₂ sequestration and natural gas storage operations, non-conventional hydrocarbon production, etc.). These activities can alter the stress field of the shallow Earth's crust by pore pressure changes or volume and/or mass changes inducing or triggering seismicity (Ellsworth, 2013; Giardini, 2009; Mignan, 2016). Such earthquakes are a nuisance or even a danger to the local population and can strongly undermine the societal acceptance of a project (Trutnevyete and Wiemer, 2017; Grigoli et al., 2017; Hirschberg et al., 2015). The recent Pohang earthquake with a magnitude of 5.5 (Grigoli et al., 2018; Kim et al., 2018) is an extreme example of a triggered earthquake related to geothermal activities that had a combined economic impact of USD 300 million as well as more than 135 injuries (Lee et al., 2019). In particular, fluid injection or extraction in tectonically active zones carries a risk of inducing a seismic event of a significant magnitude (Grigoli et al., 2017), and deep geothermal projects are a primary example. Another source of concern stems from the fact that deep geothermal projects in Europe – the Geldinganes stimulation is no exception – are often located close to consumers and thus in densely urbanized areas with historical and vulnerable buildings and infrastructures. In these contexts, the problem of assessing and managing induced seismicity is critical (Bommer et al., 2015; Giardini, 2009; Majer et al., 2007, 2012; Mignan et al., 2015, 2019a, b, Trutnevyete and Wiemer, 2017; van Elk et al., 2017; Walters et al., 2015). It is also a well-known fact that societal acceptance of induced seismicity has substantially decreased in some countries in the past decade, a result of failures discussed widely in the media and an overall change in risk perception.

Despite the large body of research conducted over the past decades by numerous research groups, the physical, chemical, and hydro-mechanical mechanisms governing induced seismicity are far from being fully understood, posing clear limits to the risk assessment and management strategies (Yeck et al., 2017; Trutnevyete and Wiemer, 2017; Grigoli et al., 2017; Mignan et al., 2019a, b). The limitations to forecasting induced seismicity are, on the one hand, the non-uniqueness on the physical framework for modeling and, on the other hand, even more importantly, the large uncertainties on the boundary conditions needed for forecasting (e.g., where are faults, what are their sizes and stress state, what is the permeability distribution of the reservoir, etc.). It follows that any risk assessment and management strategy must capture the existing uncertainties and lack of knowledge, requiring a probabilistic approach that explicitly considers both epistemic uncertainties and aleatory variability. It also implies that in order to reduce uncertainties, the risk assessment should be updated as soon as new data become available during the drilling and stimulation phase.

Despite these challenges, geothermal energy is a highly important renewable energy resource with a low carbon footprint. It has been successfully operated in many areas for decades, and Iceland is a prime example for economically successful and widely accepted use of deep geothermal energy. Several, past geothermal projects have been successfully managed with classical traffic light approaches (Majer et al., 2007; Bommer et al., 2006; Kwiatek et al., 2019) and simplified risk assessments. However, classic traffic light systems (TLSs) are simple heuristic methods, often based on expert opinions, and their likelihood of success in mitigating seismic risk is not yet clear (Baisch et al., 2019). In fact, there are several notable cases (e.g., Basel – Majer et al., 2007, and Mignan et al., 2015; Pohang – Grigoli et al., 2018; to name a few) where classical TLSs have not been successful. In most of these cases, the main event that occurred after the project was terminated, i.e., after the mitigation strategy ceased its effects. As a consequence, we consider it important for the future development of geothermal energy near urbanized areas to move beyond the existing state of the technology and develop and implement a robust, quantitative, and coherent risk management framework during all stages of a project, including post-injection.

Within this context, this study represents, to the best of our knowledge, the first publicly available probabilistic seismic risk study prior to a deep geothermal project in Iceland (and one of the very few worldwide). We attempted to combine all available risk-related information on the upcoming stimulation of the RV-43 well on Geldinganes into one quantitative and risk-based assessment. This a priori study, then, represents the basis for risk updating once the project has started and in situ real-time data become available. This procedure ideally enables a dynamic risk management solution that will also help to ensure public acceptance and thus contribute to the continued successful use of deep geothermal energy resources in Iceland and beyond. In details, the key objectives of this a priori study are as follows.

• Interdisciplinarity-based risk. Hydraulic reservoir modeling, empirical data of past sequences, expert knowledge, ground motion prediction equations, and first-order exposure and vulnerability information are integrated into one quantitative risk assessment.

• State of knowledge. Methodologies are used that are well aligned with the good-practice recommendation of the DESTRESS project (Grigoli et al., 2017; Pittore et al., 2018), with Swiss good practice recommendations (Trutnevyete and Wiemer, 2017), and with the recommendations of the international expert committee investigating the Pohang earthquake (Lee et al., 2019).

• Explicit uncertainty treatment. The uncertainties in knowledge and the variability in the data are consistently treated via use of a logic tree approach. This reflects the current state of practice in probabilistic seismic hazard and risk assessment for natural earthquakes.

• Transparency and reproducibility. The study documents all decisions taken in a transparent and reproducible way. All stakeholders in risk governance thus have access to the same level of information as a baseline and ideally a common understanding of the project's risks.

• “Updatable”. Most important, the a priori risk assessment can be updated in a consistent way as soon as new data arrive. Because the initial uncertainties are very large, updating it with in situ information is a must and should be done in a manner which is fully compatible with the initial risk assessment. The a priori risk assessment presented here is thus also a first and critical step toward risk management.

• Limitations. This study is conducted without local data and therefore uses only “off-the-shelf” models and data from different projects. It is therefore essential to clearly outline the assumptions and simplifications so that in the “update and review” phase these limitations can be addressed and the overall approach improved through more sophisticated modeling.

We fulfill these objectives by structuring the paper as follows. Sect. 2 describes the site, the geological conditions, and the planned field operations. Sect. 3 introduces the probabilistic fluid-induced seismic hazard assessment and Sect. 4 the probabilistic fluid-induced seismic risk assessment. Sect. 5 discusses the hazard and risk results as well as known limitations.

2 Site description, geological conditions, and planned operations

2.1 Site description

Well RV-43 is located on the peninsula of Geldinganes in the northeastern part of the city of Reykjavik (Fig. 1). OR is the main supplier of heat in Reykjavik and has drilled several wells on Geldinganes. It aims to produce hot water from RV-43 to be directly utilized for heating purposes and to meet the increasing energy needs of Reykjavik.

RV-43 was drilled in 2001; it is 1832 m long, where the last 1130 m is uncased (8 1∕2 in. open hole). The well is deviated towards N20^∘ E (on average), and it reaches ∼1550 m true vertical depth (TVD). The well is oriented towards the northeast of Geldinganes, an area with exceptionally high geothermal gradients that is closer than the rest of Geldinganes' wells to the extinct central volcanic system north of Reykjavik and to a possible fault zone (Steingrímsson et al., 2001). Both temperature logs and magnetic measurements support this hypothesis. Except for minor losses close to the bottom of the well, no mud losses were observed during drilling of the open-hole section of the well. The location of the well RV-43 is shown in Fig. 1.

The first and only stimulation of RV-43 took place in 2001 after its drilling. Water of pressure up to 10 MPa was injected along the open-cased segment of the well, and the total injected volume was not documented. However, this can be inferred from the original drilling report, which states that at least 1900 m³ was injected and no seismicity observed (suggesting a maximum magnitude threshold M<2, which was the minimum detectable magnitude). After the stimulation, the well had an injectivity index less than $\mathrm{6}\times {\mathrm{10}}^{-\mathrm{9}}$ m³ Pa⁻¹ s⁻¹ for the maximum injection's pressure, which is at best half of the required value for commercial exploitation.

2.2 State of stress and structural geology

A first estimate of the state of stress at Geldinganes has been inferred from a global Icelandic stress survey conducted by Ziegler et al. (2016) and by Heidbach et al. (2016). They suggest a potential orientation for ${\mathit{\sigma }}_{H_{\max }}$ of 340–40^∘ NW–SE, based on 12 geological indicators in a 10 km region around the site. The magnitude of the stress at depth could be extrapolated from shallow hydrofracturing stress measurements. Such tests were conducted in two boreholes (H32 and H18) near Reykjavik on the flank of the Reykjanes–Langjökull continuation of the Mid-Atlantic Ridge (Haimson and Voigt, 1976; Haimson, 1978).

Four tests were conducted in the borehole H32 between 200 and 375 m depth in jointed basalt. The minimum compressive stress, ${\mathit{\sigma }}_{H_{\min }}$, was found to be horizontal in the range of 4 to 6 MPa, while for the maximum horizontal stress, ${\mathit{\sigma }}_{H_{\max }}$, it was approximated as varying between 5 and 10 MPa for the four tests. The direction of ${\mathit{\sigma }}_{H_{\max }}$ was calculated based on three hydro-fractures with an orientation of N25^∘ W ± 5^∘. The vertical stress, σ_V, was calculated based on a gradient of 27 MPa km⁻¹. These values, if extrapolated to a depth of 1.5 km, suggest a normal stress regime. In the borehole H18 only three tests were performed due to extensive jointing. While the test at 180 m was conducted in basalt, the lower tests at 290 and 324 m were in an intrusive dolerite. Also, in this case, the minimum principal stress was found to be horizontal (${\mathit{\sigma }}_{H_{\min }}$), increasing with depth from 4 to 8 MPa. For the maximum horizontal stress (${\mathit{\sigma }}_{H_{\max }}$) it was observed in a range from 12 to 16 MPa, while the vertical stress ranged from 5 to 9 MPa. For these tests, the hydro-fractures suggest contradictory directions and hence two possible orientations for ${\mathit{\sigma }}_{H_{\max }}$: N20^∘ E for the 180 m test and N45^∘ W for the 290 m test. Extrapolation of the results at 1.5 km depth suggests in this case a strike–slip to the reverse regime.

Combining the results of both boreholes, a linear approximation of the data between 200 and 350 m depth gives ${\mathit{\sigma }}_{H_{\min }}=\mathrm{21}$ MPa km⁻¹, ${\mathit{\sigma }}_{H_{\max }}=\mathrm{3}$ MPa + 30 MPa km⁻¹, and σ_V=27 MPa km⁻¹. As reported by Haimson and Voigt (1976), the measured stress orientation (H32) has no obvious relationship to the NE strike of individual rift zone fissures and faults, inferred WNW direction of lithospheric plate motion, or axial rift zone earthquake focal solutions which indicate NW trending. The measured stresses could be related to (i) a hot spot, (ii) local phenomena involving the extinct NNW-trending Kjalarnes central volcano, or (iii) ground distortion due to fluid withdrawal from the Laugarness hydrothermal system. Finally, the two stress measurements in dolerite in borehole H18 could be interpreted as high-stress layers. By excluding these two measurements, the linear approximation gives ${\mathit{\sigma }}_{H_{\min }}=\mathrm{2}$ MPa + 10 MPa km⁻¹, ${\mathit{\sigma }}_{H_{\max }}=\mathrm{3}$ MPa + 13 MPa km⁻¹, and σ_V=27 MPa km⁻¹, leading to normal conditions at 1.5 km depth (Hofmann et al., 2020).

2.3 Planned activity

The re-stimulation of well RV-43 is foreseen by the end of October 2019. The re-stimulation is based on a three-staged cyclic pulse stimulation that will last for circa 12 d (4 injection days per stage). In particular, it is expected to enhance productivity in three pre-existing fracture zones penetrated by RV-43 and isolated with straddle packers. Specifically (i) the first zone is located at 1700–1750 m measured depth (MD) that corresponds to 1467–1507 m in TVD, where the basalt intersects with the gabbro and mud losses had been observed, (ii) the second zone is located at 1300–1350 m MD (1150–1189 m TVD), and (iii) the third zone is located at depth 1100–1150 m MD (1001–1032 m TVD). Each of the zones is stimulated with a cyclic injection scheme (“cyclic” stimulation), which repeats every 24 h and includes pressurizing RV-43 with pulses of frequency of 1∕60 Hz (“pulse” stimulation) and continuous injection phases. This is illustrated in Fig. 2.

Figure 2 Example of the main stimulation for one stage. After stimulation, flowback is performed.

The application of short-term cycles is based on the concept of fatigue-hydraulic fracturing, introduced by Zang et al. (2013, 2017, 2019). In practice, pressure pulses are expected to weaken the rock (“fatigue”) by inducing microcracks before macroscopic failure. This mechanism has three major intended benefits. First, the stimulated reservoir volume is increased due to more complex fracture growth, and a larger and denser fracture network provides a larger heat exchanger area. Second, the breakdown pressure is reduced, and, therefore, lower injection pressures are required to stimulate the target formation, hence reducing the potential for slip on faults and, thus, the likelihood of induced seismic events. Third, the magnitude of the largest induced seismic events is potentially limited. The stimulation of each stage is expected to last a maximum of 4 d with the following schedule, including pre- and poststimulation operations:

• 1∕2 d to install the packer,

• 1∕2 d to perform injection tests with a stepwise flow rate increase,

• 4 d for the main stimulation (with stepwise flow rate increase, and repeating phases of cyclic injection, cyclic-pulse stimulation, and continuous injection),

• 1∕2 d for performing flowback, where withdrawn water goes to the sea,

• 1∕2 d for removing the packer and redressing it for the following stage.

Injected water is not expected to exceed rates of 60 L s⁻¹ or overpressures of 20 MPa at any time during the stimulation due to restrictions by the equipment. No threshold value for injectivity has been reported for stopping the stimulation, and stimulations are expected to continue either until the end of the planned injection or until a TLS forces the termination (e.g., Mignan et.al., 2017). The well is stimulated sequentially from bottom to top.

An exemplary main stimulation for each stage is plotted in Fig. 2. During the first day, step rate injection tests are performed for estimating the pressure at which fractures open and to observe the seismic response to increasing flow rates. Based on this, the flow rates of the following phases are determined in order to reach sufficient pressures for stimulation of the target interval. The main part of the stimulation consists of cyclic injection (four cycles of 1 h high-rate injection and 1 h low-rate injection), cyclic pulse injection (four cycles of 1 h high-rate injection with pressure pulses and 1 h low-rate injection), and 8 h of continuous injection. The volume injected in each of the phases is planned to be approximately equal. The flow rates depend on the fracture opening pressure. This procedure is repeated up to three times before the flow rates are reduced slowly and stepwise at the end of the treatment.

2.4 Mitigation strategy

In the presence of fluid-induced seismic risk, it is paramount to efficiently monitor the induced seismicity and define a risk mitigation strategy. In the Geldinganes area, a dedicated microseismic network has been recently installed. The seismic monitoring infrastructure, completed in August 2019, consists of 13 seismic stations, one seismic array of 7 seismic stations, and one deep borehole array of 17 geophones (Fig. 1). The stations send data in real time to the Iceland GeoSurvey (ISOR), which streams them both to ETH Zurich and GFZ Potsdam. Real-time seismic data analysis will be performed using the software package Seiscomp3. Induced seismicity monitoring and risk mitigation operations at Geldinganes are conducted by a team of experienced professionals, including seismologists, field operation managers, reservoir engineers, and an internal expert panel (who will support decision-making during critical situations). The adopted protocol for the Geldinganes TLS is based on a five-step action plan that governs the fluid injection operations illustrated in Fig. 3 and summarized below.

Figure 3 The classic traffic light scheme adopted in the Geldinganes project.

In the case of induced seismic events above a certain threshold, we require a specific action plan. In particular, we first subdivide the region surrounding the industrial site into an internal and external domain.

• Internal domain. This defines the volume surrounding the industrial operations where seismicity will be monitored and analyzed with maximum sensitivity.

• External domain. This is a wider volume surrounding the internal domain, where the occurrence of seismicity may still be associated with the industrial operations.

For the Geldinganes site, we set these domains as cylinder-shaped volumes, with radii from the injection well of 2.5 and 5.0 km for the internal and external domain, respectively (Fig. 1). The range of depth considered for both domains is between 0 and 10 km. These values have been defined by considering other induced monitoring projects and considering the expected uncertainties for the automated locations. The magnitude of completeness of the Geldinganes network, evaluated using the BMC method (Mignan et al., 2011; Panzera et al., 2017), is ∼0.3 and ∼0.0 for the external and internal domain, respectively. Seismic events with M_L>0.0 and occurring within the internal domain should be seen clearly by almost all the stations within this domain. Therefore, all the seismic events above this magnitude threshold will be manually analyzed.

For the external domain, we will manually refine the automated solutions for seismic events with M_L>0.5. Since automatic magnitudes of small events might be overestimated, these thresholds need to be revised by considering the observed seismicity data collected during the early stage of stimulation operations. Then, the analyzed data are used to update the risk study and to assess the performance of the monitoring network. These analyses are performed at the early stage of the cyclic stimulation, and injection is increased carefully until at least a few events are detected and located.

3 Probabilistic fluid-induced seismic hazard assessment

Probabilistic risk assessment is emerging as the standard approach to manage and mitigate induced seismicity linked to fluid injections in the underground (Mignan et al., 2015, 2017, 2019a, b; Bommer et al., 2015; Grigoli et al., 2017; Lee et al., 2019). The need for a probabilistic risk-based approach is motivated by the stochastic nature of earthquakes; the many uncertainties associated with the process of inducing seismicity; and the needs of regulators, insurance, and public engagement (Mignan et al., 2019a, b). Both hazard and risk approaches follow standards proposed, among others, by the Swiss Seismological Service (Wiemer et al., 2017) and related references (Broccardo et al., 2017a; Mignan et al., 2015, 2017, 2019a, b), which are based on a combination of probabilistic seismic hazard analysis (PSHA) and the PEER-PBEE framework (Cornell, 1968; Cornell and Krawinkler, 2000).

PSHA is assessed as the probability of exceeding a given intensity at a given distance R from the injection site, based on the number of events above a given minimum magnitude m₀, the frequency distribution of the magnitude (namely the truncated Gutenberg–Richter distribution), and an empirical ground shaking attenuation function. The latter can be an intensity prediction equation (IPE) based on felt intensity or a ground motion prediction equation (GMPE) based on peak ground acceleration (PGA), spectral acceleration (SA), or peak ground velocity (PGV). Commonly, within this probabilistic framework, there are two main elements to be defined: (i) the probabilistic characterization of the seismogenic source model(s) and (ii) the ground motion characteristic model(s) (describing the expected ground vibration given the occurrence of an earthquake). The first gives the temporal and spatial forecast of the earthquake ruptures, while the second is characterized by GMPEs to link the earthquake rupture with the expected ground shaking at the site of interest.

The output of PSHA analysis is the rate of exceedance or hazard curves (probability of exceedance for a given period of time) of a given ground shaking intensity measure (IM) type. A single curve (for a given set of parameters) represents the aleatory (irreducible) variability within the defined model. To include also the epistemic uncertainties, given the alternative possible models, a logic tree structure with weighted branches (indicating the belief in a given model) is defined (e.g., Mignan et al., 2015). Figure 4 shows the proposed logic tree adopted for this a priori risk analysis. The first level of the logic tree defines the seismogenic source models, the second level the upper bound of the Gutenberg–Richter distribution, and the third level the GMPEs and the ground motion intensity conversion equations (GMICEs). In the following, we report a detailed discussion for each level of the logic tree.

Figure 4 Logic tree for the PSHA analysis. The weight is reported in the grey boxes (e.g., in the SM1 model each a_fb−b combination has uniform weight of 1∕13, each GMPE has a uniform weight of 1∕7, and each GMICE has uniform weight of 1∕2).

3.1 Seismogenic source models

In this analysis, we assume that induced seismicity nucleates and eventually extends in the proximity of the injection point. Therefore, a point source located at the coordinates of the injection point is used as the unique seismogenic source model for the investigation. This implicitly excludes any geometrical uncertainty on the location of the hypocenter. Forecasting the number of events that will occur in a reservoir stimulation is difficult because (as previously stated) the stressing conditions and location of faults near the injection point are unknown. Empirical data from similar sites can be used as a first-order proxy, but in the case of Geldinganes, only limited experience exists. In light of these limitations, we argue that the spatial variability in the seismicity is well constrained by a simple seismogenic source that can be updated for real-time application.

The number and size of earthquakes in PSHA analysis are based on three parameters that describe the local seismic activity rate, the event size distribution, and the largest event size (Cornell, 1968). These parameters are typically constrained based on observed seismicity, with the activity rate broadly scaling with the seismotectonic strain input. For induced seismicity, the seismogenic models must also describe the local seismic productivity that is (in this case) linked to the injection profile. Such seismicity rates are unknown, albeit the hydraulic energy input might be estimated beforehand. Moreover, the link between induced seismicity and stress release is a key factor to be considered in the analysis. The fraction of seismic to hydraulic energy may thus vary from zero (no events observed) to well above 1 (sometimes referred to as “triggered” events that release mostly pre-accrued tectonic stresses; e.g., Pohang – Grigoli et al., 2018). We will consider two simple seismogenic source models to analyze this uncertainty in energy release and have a first-order forecast of the underground response to injection.

• Model SM1. This is a seismogenic source model that assumes the underground feedback is site-specific constant, with all parameters purely data-driven (Dinske and Shapiro, 2013; Mignan et al., 2017; Broccardo et al., 2017a).

• Model SM2. This is a seismogenic source model that simulates the fluid and overpressure propagation for the planned injection protocol based on one-dimensional diffusion and stochastically distributed seeds. Model SM2 will also explicitly use the observation of no seismicity at M≥2 during the first stimulation in the year 2001 as a constraint. The synthetic catalogs are then converted onto the same underground feedback site-specific parameters of model SM1 (Karvounis et al., 2014; Karvounis and Jenny, 2016).

These two models capture, to a first-order, the epistemic uncertainty in forecasting seismicity, since they express alternative approaches to forecasting (purely empirical and partially physics-based). Both models are equally weighted to estimate the ground shaking estimation at the site of interest.

3.1.1 Seismogenic source model SM1

The seismogenic source model SM1 assumes that the “seismic underground feedback” per volume affected by significant pore-pressure change is a site-specific (and generally unknown a priori) constant. This constant can vary by several orders of magnitude between sites. Because the volume affected scales with the volume of fluid injected (and in theory to the pressure applied; Mignan, 2016; Langenbruch et al., 2018), this implies a relation between the expected number of earthquakes E[N] and the volume injected V, as

$$\begin{array}{}\text{(1)}& \begin{array}{rl}E& \left[N\right(t);M>m]\\ & =\left\{\begin{array}{ll}{\mathrm{10}}^{a_{\mathrm{fb}}-bm}V\left(t\right),& t\le T_{\mathrm{in}},\\ {\mathrm{10}}^{a_{\mathrm{fb}}-bm}\mathit{\tau }\exp \left(-\frac{t-T_{\mathrm{in}}}{\mathit{\tau }}\right)\dot{V}\left(T_{\mathrm{in}}\right),& t>T_{\mathrm{in}},\end{array}\right.\end{array}\end{array}$$

where a_fb is the underground feedback parameter (i.e., the overall activity for a given volume V, which is better known as the seismogenic index Σ; e.g., Dinske and Shapiro, 2013), b is the slope of the Gutenberg–Richter distribution, T_in is the injection duration, and τ is the mean relaxation time of a diffusive process. This linear relation (during the injection phase), first described by Shapiro's group, is broadly accepted in the technical community as a first-order model (e.g., Dinske and Shapiro, 2013; van der Elst et al., 2016; Mignan, 2016; Broccardo et al., 2017a). We should mention that we use the generic term a_fb instead of Σ to remain agnostic as to the physical origin of this linear relationship. The seismogenic index infers a poroelastic origin (e.g., Shapiro and Dinske, 2009), although other drivers, such as overpressure field geometry, can also explain the linearity observed between V and N (Mignan, 2016). The post-injection phase has been added by Mignan et al. (2017) to account for the decrease of the rate of seismicity after the injection has been terminated (trailing effect) (Mignan et al., 2015). This model has been verified for a number of fluid injection experiments in terms of flow rate $\dot{V}$ versus induced seismicity rate λ(t, M>m) (Mignan et al., 2017). This allows a refined analysis by defining the rate function

$$\begin{array}{}\text{(2)}& \mathit{\lambda }(t,M>m)=\left\{\begin{array}{lcl}{\mathrm{10}}^{a_{\mathrm{fb}}-bm}\dot{V}\left(t\right),& \mathrm{for}& t\le T_{\mathrm{in}},\\ {\mathrm{10}}^{a_{\mathrm{fb}}-bm}\dot{V}\left(t\right)\exp \left(-\frac{t}{\mathit{\tau }}\right),& \mathrm{for}& t>T_{\mathrm{in}},\end{array}\right.\end{array}$$

which allows the definition of a non-homogeneous Poisson process (NHPP). Note that this model only applies to the stimulation phase in which the fluids injected are not supposed to be produced back, hence creating an overpressure field at depth z. In practice, after each stimulation stage parts of the injected fluid will be produced back by natural bleed off (without pumping) directly after each stage and by airlift testing after the end of the last stage. Note also that $E\left[N\right(t),M>m]=\underset{\mathrm{0}}{\overset{t}{\int }}\mathit{\lambda }(t,M>m)\mathrm{d}t$ and the expected total number of fluid-induced earthquakes is $E\left[N\right(\mathrm{\infty });M>m]={\mathrm{10}}^{a_{\mathrm{fb}}-bm}\left(V\right(T_{\mathrm{in}})+\mathit{\tau }\dot{V}(T_{\mathrm{in}}\left)\right)$.

While the parameters a_fb, b, and τ can be estimated during the stimulation (Mignan et al., 2017; Broccardo et al., 2017a), a priori knowledge on those parameters is limited and the range of possible values wide. Given this state of “uniform” uncertainty, we assigned equal weights to all the possible [a_fb, b] combinations. We list the a_fb, b, parameter estimates for different sites in Table 1, which will be used as input for the a priori risk study. Uncertainties are likely to significantly reduce once seismic data are obtained by monitoring during the stimulation. Note that due to the correlation of a_fb and b (Broccardo et al., 2017a), pairs of (a_fb, b) values from different sites need to be maintained. In Mignan et al. (2017), the mean relaxation time has been observed to widely vary between injection sites with $\mathrm{0.2}<\mathit{\tau }<\mathrm{15}$ d.

Table 1 Underground seismic feedback to deep fluid injection.

¹ ISO code. ² Referred to as seismogenic index in Dinske and Shapiro (2013).

In order to apply a classical PSHA analysis, we transform the NHPP into a homogeneous Poisson process (HPP), using the equivalent rate ${\mathrm{\Lambda }}_{M>\mathrm{2}}=E\left[N\right(T),M>\mathrm{2}]=\underset{\mathrm{0}}{\overset{T}{\int }}\mathit{\lambda }\left(t\right)\mathrm{d}t$, for a unit of time which corresponds to the total project period (including the post-injection phase); i.e., $T=T_{\mathrm{in}}+T_{p_{\mathrm{in}}}$, where $T_{p_{\mathrm{in}}}$ is the post-injection time. We selected M>2 because we assume that lower magnitudes will not have the potential to trigger any damage. By doing so, $P(M>\mathrm{2};T)=\mathrm{1}-\exp (-{\mathrm{\Lambda }}_{M>\mathrm{2}})$. Table 1 and Fig. 5 (red dots) report the equivalent rate Λ_M>2 for each project, for a target injected volume of circa V=18 000 m³ (estimated from ca. 6000 m³ injection per stimulation, multiplied by 3 stimulations; Fig. 2). At present, without any pre-stimulation phase it is not possible to infer where the Geldinganes project is placed in this domain; however, what is known is that 5000 m³ of water was injected and no seismicity observed.

3.1.2 Seismogenic source model SM2

With model SM2, a first-order physical process is included into the forecasting. This is done by modeling pressure diffusion through a fractured media containing randomly distributed earthquake faults (so-called “seeds”). The pressure propagation can be adopted based on the reservoir properties, limited to the available information. Then, the density of these seeds and their size distribution are treated as free site-specific parameters that (again) are unknown a priori. These models are commonly referred to as “hybrid” models (Gischig and Wiemer, 2013; Goertz-Allman and Wiemer, 2013), as they combine deterministic and stochastic modeling. Specifically, the adaptive hierarchical fracture representation (a-HFR) is employed both for modeling flow in a fracture network with dynamically changing permeability (Karvounis and Jenny, 2016) and for simulating the source times of randomly pre-sampled scenarios of hydro-shearing events at certain hypocenters (Karvounis et al., 2014). This hybrid model is chosen here, as it can integrate several of the field observations, returns forecasts both of the spatial distribution of seismicity and of its focal planes, and can forecast reservoir properties like the expected well's injectivity at the end of the injection.

The required inputs for the proposed hybrid model are the initial hydraulic properties, the planned activities, a first-order knowledge of the stress conditions in the proximity of the well, and the orientations of pre-existing fractures. Here, the extrapolated stress measurements described in Sect. 2.2 are employed without excluding any of the measured stresses; i.e., the vertical direction is a principal direction and ${\mathit{\sigma }}_{H_{\min }}=\mathrm{21}$ MPa km⁻¹, ${\mathit{\sigma }}_{H_{\max }}=\mathrm{3}$ MPa + 30 MPa km⁻¹, and σ_V=27 MPa km⁻¹. The planned activities are described in Sect. 2.3, the initial transmissibility is in agreement with the currently expected injectivity of $\mathrm{6}\times {\mathrm{10}}^{-\mathrm{9}}$ m³∕(Pa s⁻¹) along the whole open segment of RV-43, and the compressibility is inferred from the nearby well HS-44, since there are no reported measurements of the characteristic time at RV-43. Moreover, in this a priori analysis, all surface orientations are considered equally probable. Observe, however, that two distinct sub-vertical fault sets seem to prevail at the mainland surrounding the bay above the injection stages (Hofman et al., 2020).

We assume a constant value for those parameters for which the rate of seismicity is less sensitive. These are as follows: the friction and the cohesion of fractures for the Mohr–Coulomb failure criterion which are fixed equal to 0.6 and 0, respectively, and the mechanical aperture of fractures after they have slipped is 1 mm. The remaining parameters (those for which seismic rate is more sensitive) are the spacing between pre-existing fractures, the b value of the Gutenberg–Richter during injection, and the permeability of a slipped surface. A range of possible values is considered for each of the latter parameters; i.e., the spacing between pre-existing fractures (1, 2, 4, 8, 16 m), the b value of the Gutenberg–Richter law (0.75, 1.0, 1.25, 1.5, 1.75, 2.0), and the hydraulic aperture between the two parallel fracture surfaces (200, 500 µm), out of which the equivalent permeability can be estimated. A synthetic catalog of rate Λ_M>2 is created for each possible combination of the above values, and an equivalent Λ_M>2 is computed. Next, following the same logic of Sect 3.1.1, we assign equal weights to each a_fb−b combination. The resulting scenarios are shown in Fig. 5 as blue dots in a_fb−b space. Note that all scenarios are above the dashed line that indicates the limit posed by the no-seismicity observation during the initial stimulation.

Figure 5 Distribution of a_fb−b values for the synthetic catalog together with the dataset of Table 1. In the planned injection profile (see Fig. 2), the flow rate decreases progressively back to zero, meaning that this simple model cannot strictly be applied. As approximation, we use max(ΔV)=1728 m³ d⁻¹ instead of ΔV_shut-in. A direct comparison can be made between the volume injected V=18 000 m³ and the equivalent τΔV=2880 m³ for τ=1 d and 28 800 m³ for τ=10 d. The dashed line represents the upper limit of no expected seismicity of M>2.

3.2 Upper bound for the Gutenberg–Richter distribution vs. maximum observed magnitude distribution

The frequency magnitude distribution of natural and induced earthquakes follows (to a first order) the classical Gutenberg–Richter distribution. This distribution is truncated at an upper end for energy conservation but also because existing faults and fault systems have a maximum size (e.g., Bachmann et al., 2011; Mena et al., 2013; Mignan et al., 2015; Baker and Gupta, 2016). Empirical observations will only poorly constrain the largest possible earthquake, since it is by definition an exceptionally rare and extreme event. One of the major sources of uncertainty is thus in PSHA related to the upper bound of the (truncated) Gutenberg–Richter distribution, here indicated as m_max.

It is generally accepted that the largest possible induced earthquake cannot be larger than the tectonically largest one. However, in induced seismicity, the tectonic environment (controlled primarily by the state of stress) at a site may be such that no tectonically pre-stressed larger ruptures exist. Under these conditions, ruptures will run out of energy once they leave the volume brought into a critical state for failure by the injection – e.g., because of the effect of pore pressure on the Coulomb failure criteria. In these conditions, run-away ruptures cannot occur even if a natural fault exists (in other words, triggered events cannot happen), and the largest magnitude size, as a consequence, is limited by the volume or area affected by overpressure (which again scales with the volume of fluid injected and the hydraulic properties of the subsurface). In such a situation, m_max can locally be substantially smaller than the regional tectonic one. This is common in “fracking” operations in tight shales. McGarr (1976, 2014), formalized this volume limit as $m_{\max ,\mathrm{McGarr}}=\mathrm{2}/\mathrm{3}{\log }_{\mathrm{10}}\left(GV\right)-\mathrm{10.7}+\mathrm{14}/\mathrm{3}$, where $G=\mathrm{3}\times {\mathrm{10}}^{\mathrm{10}}$ Pa is the modulus of rigidity. McGarr has shown that this relationship is consistent with the data from a compilation of injections. However, a number of researchers (Gischig et al., 2014; van der Elst et al., 2016; Mignan et al., 2019b) have pointed out that outliers exist (e.g., Pohang, South Korea – Grigoli et al., 2018, and Kim et al., 2018; St. Gallen, Switzerland – Diehl et al., 2017) and that the McGarr limit is best explained as a purely statistical relationship based on simple extreme-value-theory principles (Embrechts, 2013).

The McGarr limit has been used (and in some cases one might argue misused) in numerous induced seismicity hazard assessments (van der Elst et al., 2016). For V=18 000 m³, we would for example obtain $m_{\max ,\mathrm{McGarr}}=\mathrm{3.79}$. Based on the recent statistical tests of van der Elst et al. (2016), and the occurrence of the 2017 Pohang earthquake above the expected limit (Grigoli et al., 2018; Kim et al., 2018), a fixed McGarr limit now appears questionable to many seismologists. Therefore, since the complete information about the number, location, size, and stressing condition of faults in the Geldinganes area is not available (in particular before the stimulation phase), it is appropriate to consider m_max to be the regional tectonic m_max=7 (Kowsari et al., 2019). This estimate could be reduced at a later stage if local fault information were found to provide better constraints. In fluid-induced seismicity, when m_max is related to the tectonically largest event, it is not a critical choice (Gupta and Baker, 2017). This is because the rate of occurrence is typically significantly low compared to the typical return periods of interest. It follows that both hazard and risk generally are dominated by the more frequent, moderately sized events.

In addition, in this study, we determine the probability distribution of the maximum observed magnitude, M_T¹, at a fluid injection sites for the total time of observation T (Holschneider et al., 2011). This is fundamentally different from the upper bound m_max of the Gutenberg–Richter distribution, which is merely a deterministic upper limit fixed by physical constrains. The probability distribution of the maximum magnitude, M_T=max[M₁, …, $M_{t_i}$, …, M_T], can be easily derived considering that the magnitude events are statistically independent. It follows that $F_{M_T}(m\mathrm{|}N=n)=F_M(m{)}^n$, where F_M(m) is the classical Gutenberg–Richter cumulative probability density function, and $f_{M_T}(m\mathrm{|}N=n)=n{F}_M\left(m{)}^{n-\mathrm{1}}{f}_M\right(m)$. Since the number of events is a random variable itself, then $F_{M_T}\left(m\right)=\sum _n{F}_{M_T}(m\mathrm{|}N=n)P(N=n\mathrm{|}\mathrm{\Lambda }(T\left)\right)$, $f_{M_{\max }}(m\mathrm{|}N=n)=\sum _{n}n{F}_M\left(m{)}^{n-\mathrm{1}}{f}_M\right(m\left)P\right(N=n\mathrm{|}\mathrm{\Lambda }\left(T\right))$, where P(N=n|Λ(T)) is the classical Poisson discrete distribution. Figure 6a and b show the equivalent rate of seismicity, $\mathrm{\Lambda }\left(T\right)(\mathrm{1}-F_{M_{\max }}(m\left)\right)$ (i.e., a weighted complementary CDF) for each of the projects reported in Table 1 (SM1 model) and for each of the synthetic catalogs (model SM2). We can observe a large scatter of the rate of seismicity, reflecting the large uncertainty exiting prior to the project.

Figure 6 Envelope probabilistic density distribution of the rate model and maximum observed magnitude M_max (a, c) based on Table 1 and (b, d) based on synthetic catalog (S2 source model). (e) Expected magnitude per volume injected based on Table 1.

Together with the distribution of M_T for each a_fb−b couple, we report the envelope distribution computed as the mean value over all the branches of the logic tree (Fig. 4). Figure 6c shows the envelope distribution. Observe that given the sparse dataset (Table 1), this distribution is (inevitably) multimodal. The expected E[M_T], based on this envelope distribution, is 2.25, and the 5 %–95 % interval is [0.10–4.45]. It is important to highlight that these values represent some statistics based on previous projects and not the expected values for this project. In fact, here, the envelope distribution represents a prior distribution, which must be updated during a pre-stimulation phase and during the stimulation. In the following, we also report the envelope distribution of M_T based on the synthetic catalog derived according to the SM2 source model. Figure 6d shows the envelope distribution. Different from the envelope distribution based on the SM1 source model, this distribution shows a more regular shape, since the synthetic dataset is denser and more confined. However, this prior distribution can be affected by overfitting, since it is based on stress measurements (without considering their uncertainties) that might not represent the current local condition correctly. The expected E[M_T], based on this envelope distribution, is 2.09, and the 5 %–95 % interval is [1.19–3.42]. Finally, in Fig. 6e we reported the E[M_T] and [5 %–95 %] confidence bound as a function of the injected volume.

3.3 Ground motion prediction equations and intensity measures

The relationship between the site source characteristics and given ground shaking IM types is given by seven GMPEs. Kowsari et al. (2019) provide a set of adjusted GMPEs that have been selected for this investigation. In particular, the proposed GMPEs were adjusted using newly compiled ground motion records of six strike–slip events in the South Iceland Seismic Zone (SISZ), with a range of magnitudes of M∈[5, 6.5] (M is intended as M_w, as used in Kowsari et al., 2019) and distance of R∈[0, 80] km. The intensity measures are reported in Table 2, and the value of the functional form and the coefficients can be retrieved directly from Kowsari et al. (2019). Observe that from the original list we replaced the proposed GMPE of Lin and Lee (2008) for northern Taiwan with the local GMPE (RS09), Rupakhety and Sigjörnsson (2009), which is consistent with the strike–slip nature of Icelandic earthquakes. The recalibration has been performed only for the PGA; therefore, in the following, we assume only this physical intensity measure. The selected site-to-source distance is the Joyner–Boore metric (R_JB) (i.e., the closest horizontal distance to the vertical surface projection of the fault). When the distance metric of the original GMPE is different from R_JB, the same transformations proposed in Kowsari et al. (2019) are applied. In the Supplement (Fig. S1), we show the trellis plots for the selected GMPE models.

Table 2 List of GMPEs used in this study.

It is important to highlight the limitation of these choices. First of all, the GMPEs are calibrated for natural events that are considerably larger in magnitude compared to the expected fluid-induced events (Fig. 6). Therefore, the extrapolation to lower magnitudes is biased (Bommer et al., 2007; Baltay and Hanks, 2014). This will have a more significant effect on the low-damage threshold, while the IR computations are less impacted, since they depend on larger events. Moreover, for small events in the proximity of the injection point, the ideal source-to-site distance is the hypocentral distance and not R_JB (observe that in this case R_JB converges to the epicentral distance), which neglects the hypocenter depth. As a consequence, this analysis is independent of the injection depth. Again, this limitation has an impact on the small damage threshold, since the depth of the events is expected to have a significant influence.

In this a priori assessment, we use, as a final IM, the European Macroseismic Scale (EMS98; Grünthal, 1998)². The advantage of EMS98 over instrumental intensity measures, in this phase, lies in the easier interpretability of this scale, which is based merely on shaking indicators expressed in terms of damage and nuisance to the population. Based on these considerations, the selected GMPEs are converted into expected intensity by using GMICEs for small–medium intensities. The GMICEs used in this work are introduced by Faccioli and Cauzzi (2006) and Faenza and Michelini (2010). The aleatory variability is then combined into a GMPE–GMICE model, with σ_TOT defined as ${\mathit{\sigma }}_{\mathrm{TOT}}=\sqrt{\left({\mathit{\sigma }}_{\mathrm{GMPE}}^{\mathrm{2}}\right)a^{\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{GMICE}}^{\mathrm{2}}}$, and values of the mean σ_GMPE, σ_GMICE, and a reported in the Supplement, together with the combined trellis plots (Table S1, Fig. S2).

3.4 Probabilistic hazard results (PGA, EMS98)

The hazard integral is reduced to the marginalization of the random variable magnitude, M, and the conditional random variable IM|M=m, since the site-to-source distance is fixed by the source point (which is assumed at the injection point). For a given site, then the rate of exceedance is simply reduced to $\mathrm{\Lambda }(\mathrm{im};T,b)=-\underset{m}{\int }P(\mathrm{IM}>\mathrm{im}\mathrm{|}M=m,r)d{\mathrm{\Lambda }}_{M>\mathrm{2}}(m;T,b)$, where $d{\mathrm{\Lambda }}_{M>\mathrm{2}}(m;T,b)={\mathrm{\Lambda }}_{M>\mathrm{2}}\left(T\right)F\left(m\right)$, with F(m) equal to the Gutenberg–Richter above a magnitude of 2. Given the discussion in Sect. 3.1.1, the probability of exceedance of an intensity, IM=im, for a given time period (which corresponds to the total duration of the project given the normalization introduced in Sect. 3.1.1), is given by the Poisson distribution as $P(\mathrm{IM}>\mathrm{im},t=T)=\mathrm{1}-\exp (-\mathrm{\Lambda }(\mathrm{im};T,b\left)\right)$.

Λ_M>2(T) is not known a priori (moreover, an uncertainty quantification based on local condition cannot be carried out a priori); therefore we compute the risk for each of the a_fb and b pairs of Table 1. As mentioned in Sect. 3.1.1 the scatter is very large, and this reveals the state of deep uncertainty existing prior to a pre-stimulation phase. The PSHA outputs are shown in red in Fig. 7 for both the PGA (top panels) and for the IM (bottom panels). These curves confirm the state of deep uncertainty, in particular for the location in proximity of the injection point. In fact, for a given probability of exceedance of 10⁻⁴ and distance 2–5 km from the injection point, the macroseismic intensity range between the 10 % and 90 % percentile is circa IM∈[6, 11]. In addition, we report the PSHA analysis based on the source model SM2. The outputs are shown in blue in the same panels. The epistemic uncertainties of the source model SM1 are considerably higher than the ones arising from the source model SM2. This was expected, given the inherent sparsity present in the data of Table 1. Moreover, the epistemic median of source model SM1 is higher than source model SM2. In addition to the PSHA output, in Fig. S3, we reported the hazard-based scenarios for different magnitude.

Figure 7 PSHA analysis comparison between source model SM1 (Table 1) and SM2 (synthetic catalog). Solid lines: medians. Dashed lines: 10 % and 90 % quantiles. Intensity measure: EMS98.

4 Probabilistic fluid-induced seismic risk

In seismic risk assessment, it is common to distinguish between physical and non-physical risk. Examples (and precedents) of non-physical risk include noise, vibrations felt, opposition by residents, public campaigns against the project, etc. Non-physical risk is complex and often impossible to quantify. Therefore, an effective and practical approach should focus on non-physical risk identification and mitigation rather than risk assessment (Bommer et al., 2015). Conversely, the physical risk faced by exposed communities needs a quantitative assessment. In this study, we focus only on one type of physical risk: the seismic risk.

The physical risk is commonly divided into two major categories, i.e., fatalities and/or injuries and economic losses (with both categories depending on physical damage to buildings). The a priori risk analysis for the Geldinganes project here focuses on the first risk, while the aggregate economic losses are not directly computed. Here, as a substitute for aggregate losses, we define a low-damage threshold for statistical average classes of Icelandic buildings. In particular, in this study, we select two risk measures: individual risk (IR) and damage risk (DR). IR is defined as the frequency over the time span of the project (including the post-injection phase) at which a statistically average individual is expected to experience death or a given level of injury from the realization of a given hazard (Jones, 1992; Jonkman et al., 2003; Broccardo et al., 2017b). We here define DR as the frequency over the time span of the project (including the post-injection phase) at which a statistically average building class is expected to experience light non-structural damage from the realization of a given hazard.

Since there are currently no universally used regulatory and industry approaches to manage induced seismicity of geothermal and other energy projects, we define the following safety thresholds for IR and DR. The proposed IR safety threshold is IR${}^{\mathrm{ST}}={\mathrm{10}}^{-\mathrm{6}}$. This value is more conservative compared to the typical standards for anthropogenic activities, for example in Switzerland or the Netherlands (van Elk et al., 2017). In the presence of epistemic uncertainties, the median of the IR distribution is taken as the reference metric to be compared with the selected safety standard, i.e., $q_{\mathrm{IR},\mathrm{0.5}}\le {\mathrm{IR}}^{\mathrm{ST}}$ (where q_IR,0.5 is the epistemic median of the individual risk distribution). The proposed DR threshold is DR${}^{\mathrm{ST}}={\mathrm{10}}^{-\mathrm{2}}$. As for the IR, in the presence of epistemic uncertainties, the median of the DR distribution is taken as a reference metric to be compared with the selected safety standard, i.e., $q_{\mathrm{DR},\mathrm{0.5}}\le {\mathrm{DR}}^{\mathrm{ST}}$ (where q_DR,0.5 is the epistemic median of the individual risk distribution).

The framework used for the computation of IR and DR is based on the convolution of the hazard model with the vulnerability models for the relevant building types and (only for the IR) with the consequence model. For the fragility–vulnerability model, we should base our analysis on local functions. However, at present only local fragility functions for low damage exist (Bessason and Bjarnason, 2015). Given that, we decide to use the macroseismic intensity approach for IR (Lagomarsino and Giovinazzi, 2016), while using the local fragility function for DR.

4.1 Individual risk computation

For IR, we use a vulnerability given in terms of macroseismic intensity, which follows the macroseismic approach for damage assessment (Lagomarsino and Giovinazzi, 2016), modified in Mignan et al. (2015), for the induced seismicity case. The macroseismic model defines the mean damage grade, μ_D(im), as function of a vulnerability index, V, a ductility index, Q, and a reduction factor α introduced in Mignan et al. (2015) to recalibrate low-damage states to the damage observed in the Basel 2006 sequence. The vulnerability index depends on the building class and construction specifics, and it includes probable ranges $V^{-}{V}^{+}$ as well as less probable ranges $V^{--}{V}^{++}$. Following the Icelandic exposure information described in Bessason and Bjarnason (2016), we select three building typologies, concrete, wood, and masonry, as a surrogate for pumice buildings. Moreover, Bessason and Bjarnason (2016) observed that (on average) the Icelandic buildings are stronger and more reliable than the ones based in the Mediterranean region in Europe. Based on these considerations, we select V₀ as the vulnerability index for concrete and wood and V⁻ for masonry. The choice of V⁻ for masonry is given by the observation that the fragility of this building is close to an old (before the 1980s) Icelandic reinforced concrete building. Moreover, there is no detailed information on the ductility index for the different class of building; therefore we use Q=2.3, which is the value for masonry structures and reinforced concrete structures with no seismic details. In this phase, this a practical and conservative choice, since Q=2.3 is a lower bound of the possible ranges of values for the ductility index. We report, in Table S2, the vulnerability indices together with vulnerability functions (Fig. S4) obtained by using the macroseismic model with parameters V₀=0.5 and Q=2.3.

We computed the marginal IR considering all [a_fb, b] couples in a given location (i.e., different distances), for the total duration of the project (including the post-injection phase), for both rate models SM1 and SM2, and using the HAZUS consequence model (Galanis et al., 2018; HAZUS MH MR3, 2003). The results are shown in Fig. 8 for each building class. The median and quantiles are computed considering a 50 % weight for the SM1 model and 50 % weight for SM2 model. Despite the median for each class being below the fixed threshold, i.e., $q_{\mathrm{IR},\mathrm{0.5}}\le {\mathrm{IR}}^{\mathrm{ST}}={\mathrm{10}}^{-\mathrm{6}}$, the uncertainty is very large, indicating that uncertainty quantification updates are necessary to reduce the [a_fb, b] uncertainties. Moreover, we would like to highlight that the median based merely on SM2 is considerably lower than the median based on SM1 (this is not reported in Fig. 8 for clarity). This is due to the lower variability in [a_fb, b] arising from the synthetic catalog.

Figure 8 Marginal IR for 2 and 5 km distances based on the final model (combined SM1 and SM2). The solid horizontal lines represent the weighted median values of the 1022 (13[a_fb, b]×7 GMPEs × 2 GMICE, weight 0.5, +60[a_fb, b]×7 GMPEs × 2 GMICE, weight 0.5) vertical grey lines. The dashed horizontal lines represent the 10 % and 90 % epistemic quantiles.

For DR, we use the local fragility model developed by Bessason and Bjarnason (2016). Three major categories of buildings characterize the Icelandic exposure model: reinforced concrete, timber, and hollow pumice block. Further details on the exposure model are given in Sect. S3. Within these categories, Bessason and Bjarnason (2016) define the following subcategories.

• Low-rise reinforced concrete is the first subcategory, which includes the following:

  • RCb80 – reinforced concrete structure designed before seismic code regulations (before 1980)

  • RCa80 – reinforced concrete structure designed after seismic code regulations (after 1980).

• Low-rise timber structures are the second subcategory, which includes the following:

  • Tb80 – timber structure designed before seismic code regulations

  • Ta80 – timber structure designed after seismic code regulations.

• Hollow pumice blocks (HP) are the third subcategory.

Fragility functions are provided for all these categories only for small damages (which makes the use for IR impossible). Fragility function details and damage-based scenarios for different magnitudes are reported in Sect. S3.

Finally, we computed the marginal DR considering all the [a_fb, b] couples (for both the source model SM1, with weight 50 %, and SM2, with weight 50 %) in a given location (i.e., different distances) for the total duration of the project (including the post-injection phase). The results are shown in Fig. 9 for each class of building. Again, despite the median for each class being below the fixed threshold, i.e., $q_{\mathrm{DR},\mathrm{0.5}}\le {\mathrm{DR}}^{\mathrm{ST}}={\mathrm{10}}^{-\mathrm{2}}$, there is a need for uncertainty quantification updates to reduce the [a_fb, b] uncertainties.

Figure 9 Marginal DR for the final model for 2 and 5 km distance. The solid horizontal lines represent the median values of the 511 (13[a_fb, b]×7 GMPEs, 0.5 weight, +60[a_fb, b]×7 GMPEs, 0.5 weight) vertical grey lines. The dashed horizontal lines represent the 10 % and 90 % epistemic quantiles.

4.2 Sensitivity analysis

In this section, we describe a sensitivity analysis of the epistemic uncertainties with respect to the quantities of interest (QoIs) IR and DR. Specifically, we analyze the sensitivity to the earthquake rate model and the GMPE (and GMICE), which (here) are the only source of epistemic uncertainty. The goal of this sensitivity analysis is to calculate which source of the two input uncertainties is dominant. Specifically, we performed two sensitivity analyses: one for the dataset in Table 1 and one for the synthetic catalog. This allows for better understanding of the relative contribution of the input uncertainties for each dataset.

In this study, we adopt a screening method which aims to preliminarily and qualitatively analyze the most important input parameter. In particular, we develop a modified version of the Morris method (Morris, 1991), which solves some drawbacks of the “tornado diagram” (Porter et al., 2002) used in Mignan et al. (2015). A tornado diagram is a type of sensitivity analysis based on a graphical representation of the independent contribution of each input variable to the variability in the selected QoI. Specifically, given a base model, for each considered variable, we estimate the maximum positive and negative swing of the QoI while holding all the other parameters fixed to their base value. A drawback of the method is that results are strongly dependent on the base model (i.e., it is a local sensitivity method). Therefore, we introduce a variation of the method to obtain a global sensitivity measure. The complete details of the introduced method are reported in Sect. S4, while here we discuss the general principles and the results. To obtain a global sensitivity measure, we first define a normalized local sensitivity measure of the parameter i with respect to the base model j, d_i(j) (Eq. S1). Then, we define two global sensitivity measures: the average, ${\mathit{\mu }}_{d_i}$, and the maximum, ${\stackrel{\mathrm{‾}}{d}}_i$, of d_i(j) (Eqs. S1 and S2). The sensitivity measure μ_i describes the average relative contribution of the parameter i over all possible base models j. The sensitivity measure ${\stackrel{\mathrm{‾}}{d}}_i$ describes the maximum contribution of the parameter i over all possible base models j. The two measures in this form are not normalized to 1.

Given IR and DR, Figs. 10 and 11 show the sensitivity results based on ${\mathit{\mu }}_{d_i}$ for each building class. Both measures show the same pattern. In particular, the dominating source of uncertainty is (as expected) the rate model. It is interesting to remark that the rate model contribution is more dominant for the dataset based on Table 1 than for the synthetic dataset. Moreover, this fact is consistent across different building typologies and risk metrics. This corroborates the observations that we have previously made; i.e., the uncertainties related to real data are larger than the synthetic ones (which might be affected by overfitting). The same trend is observed for ${\stackrel{\mathrm{‾}}{d}}_i$ (Figs. S8 and S9).

Figure 10 Sensitivity analysis of IR (observe that the QoI is log IR) based on the sensitivity measure ${\mathit{\mu }}_{d_i}$ for each building class.

Figure 11 Sensitivity analysis of DR (observe that the QoI is log DR) based on the sensitivity measure ${\mathit{\mu }}_{d_i}$ for each building class.

5 Discussion and conclusion

This study represents the summary of a collective effort for assessing the a priori seismic risk for the hydraulic stimulation of a geothermal well on Geldinganes, Iceland. The key findings of the assessment are shown in Figs. 8 and 9 and summarized below.

The overall risk for an individual to die in a building within a radius of 2 km around the well (Fig. 8) is assessed to be below 10⁻⁷ or at 0.1 micromort (1 micromort a is unit of risk defined as one-in-a-million chance of death). This value is within the acceptable range when compared to acceptance criteria applied in the Netherlands (or Switzerland). The reason for the acceptable risk is the overall quite limited injection volume, the fact that the initial stimulation has not produced M>2 seismicity, and the estimated low vulnerability of the building stock.

The chance of damage to buildings is around 0.1 % (Fig. 9) and therewith below the 10⁻² acceptance threshold we arbitrarily introduced for damage.

The thresholds proposed in the classical traffic light (Fig. 2) are consistent with the risk thresholds computed; it is not suggested to define more conservative TLS thresholds at this point.

The uncertainties at this stage of the project are very high, highlighting the importance of updating the risk study continuously as new data become available.

Based on the following results and the mitigation strategies summarized in the following document, we suggested proceeding with the project. However, based on the online updates of the risk model, we recommended a (possible) review of the analysis. In particular, if the median of the IR and DR grew close to the assigned limits, we prescribed a refinement of the study to address the current limitations.

5.1 Limitations of our study

Probabilistic risk assessment is in many ways a very pragmatic approach that systematically collects available information based on the current state of knowledge. It is acceptable that in many areas, the state of knowledge is limited and evolving. While we consider the current assessment to be useful and usable, there are also some limitations and areas where further improvements would be beneficial.

• Geological and seismotectonic knowledge is poorly represented. This is mostly a consequence of the fact that knowledge of the local seismotectonic condition is limited and uncertain, especially when extrapolated to the reservoir depths. The limited use of geological constraints is also a consequence of the fact that geological knowledge cannot be readily transferred into forecasting models of seismicity.

• Empirical data from similar injections in the surroundings of Geldinganes or from areas with comparable conditions are limited and mostly based on observation in the 1970s with limited seismic monitoring in place. While countless well-monitored injections have been conducted in Iceland overall, there has been less activity near Reykjavik. The initial stimulation of the Geldinganes well in 2001 produced no noticeable seismicity, which provides important constraints (Fig. 5). However, monitoring was at that time quite limited, so an event with a magnitude smaller than 2 may have been undetected, and we also need to consider that the response to the 2019 stimulation may be different.

• The seismicity forecasting models we use are simplistic in many ways, considering a limited amount of physical, hydraulic, or geological aspects. In particular, neither the SM1 nor the SM2 model explicitly consider the (re-)activation of the cracks or faults responsible for the mud losses during drilling and reported by Steingrímsson et al. (2001). We also use few models overall and do not take the risk-limiting effect of mitigation measures explicitly into account.

• Ground motion models which are specific for Iceland exist. However, they originate from few strong-motion data at short distances, from larger magnitudes, and from natural earthquakes and are therefore a limited constraint. Likewise, little is known about the site amplification at a microzonation level. We do not plan, in the advanced traffic light system (ATLS), to perform an online updating of the coefficient of the ground motion models.

• Building vulnerabilities are known at a first-order level, but no efforts have been made to verify or validate them, nor will they be updated during the ATLS implementation. No sensors in buildings are planned.

5.2 Recommendations

Below, we list a number of recommendations on risk management that were made to Reykjavik Energy before the start of stimulation operations at the Geldinganes site. These are based partially on this study but also consider experiences of past projects.

Excellent seismic monitoring and reliable near-real-time processing is a key requirement for updating the a priori risk assessment. The network installed at Geldinganes should be capable of this task; however, owing to the low seismicity in the region and short deployment time of the full network, the actual capabilities and operation procedures are untested.

The standard traffic light system operated by ISOR on behalf of OR and based on Icelandic Meteorological Office (IMO) magnitudes is critically important and the ultimate decision tool. A TLS is a simple well-proven and well-established technology; it cannot and should not at this stage be replaced with more adaptive concepts of risk assessment.

Given the uncertainties, updating this risk assessment is a key requirement. The most basic approach is to update it based on periodical reassessment of the model parameters (seismicity rates, b value, hydraulic parameters) performed offline and interactively.

A pre-stimulation test that results in a number of micro-earthquakes below a magnitude of 1, followed by a subsequent update of this risk study, would help to calibrate the seismic forecast models and to allow constraining the uncertainties. The test would also demonstrate the ability of the monitoring network to detect and locate microseismicity.

The stimulation plan consists of three distinct stages at different sections of the well. The design of this first stage is conservative, with slowly increasing flow rates and longer shut-in phases. As a consequence, the calibration of the seismic forecast models is performed along with the first stage. Then, the updated values are used as prior information for the second and third sub-stimulations.

The project is not free of seismic risk; there is a residual chance that, despite all mitigation measures applied, damaging earthquakes might occur. This report attempts to quantify this chance, and we believe it is important to openly communicate this remaining risk and the steps taken to reduce and control it to the public and authorities. This might include clarification on how potential damages would be reported, settled, and insured.

Re-activating pre-existing and tectonically pre-stressed larger fracture zones and eventually triggering a larger earthquake, like what happened in Pohang, is unlikely but still probably the most important risk for the project. The probabilistic risk approach applied here captures this chance to trigger such an event to a certain extent and in a statistical approximation. However, it may possibly underestimate the chance of such a triggered earthquake if an unknown major fault is very close to the injection site (i.e., closer than 1 km). Moreover, in this project, a fault zone may potentially be the cause for the high temperature. Observe that in Iceland the fault zones are oftentimes the targets of geothermal wells. Therefore, we suggest that the seismicity analyst team should be on the lookout for lineament, potentially indicative of a major fault zone being re-activated, and discuss it with the expert group that is accompanying the project. In particular, an in-depth analysis should be carried out after the first stage of the stimulation.

The size distribution of induced earthquakes critically determines the risk, and an unusually low b value may indicate the presence of critically stressed faults and will result in much larger probability of larger events. The reassessed b value must flow into the update of the risk assessment, but we suggest adding, as an additional safety criterion, a project halt if the b value of induced events is estimated below 0.8.

It is universally accepted that the seismicity and thus risk will decrease once the injection has been stopped. It is less clear, however, if gradual pressure reduction, shut-in, bleed off, or actively pumping out (if possible) are the best mitigation strategies. In this project, a common agreement was reached on considering bleed off to be the most adequate strategy.

Surprising developments are possible, if not likely. Therefore, we set up a small interdisciplinary expert group that can come together rapidly (e.g., virtually) if unexpected developments occur (lineaments, clusters, etc.).