Probabilistic Evaluation of Susceptibility to Fluid Injection-Induced Seismicity Based on Statistics of Fracture Criticality

Fault reactivation and associated microseismicity pose a potential threat to industrial processes involving fluid injection into the subsurface. In this research, fracture criticality, defined as the gradient of critical fluid pressure change to trigger seismicity (Δpc/h), is proposed as a novel reservoir depth-independent metric of fault slip susceptibility. Based on statistics of the fracture criticality, a probabilistic evaluation framework for susceptibility to injection-induced seismicity was developed by integrating seismic observations and hydrogeological modelling of fluid injection operations for faulted reservoirs. The proposed seismic susceptibility evaluation method considers the injection-driven fluid pressure increase, the variability of fracture criticality, and regional fracture density. Utilising this methodology, the probabilistic distribution of fracture criticality was obtained to evaluate the potential for injection-induced seismicity in both fault and off-fault zones at the Hellisheiði geothermal site, Iceland. It has been found that the fracture criticality within both fault and off-fault zones shows natural variability (mostly ranging between 0.001 and 2.0 bar/km), and that fault zones tend to be characterised by larger fracture criticality values than the off-faut zones. Fracture criticality values estimated within each zone roughly follow a Gaussian distribution. Fault zones around five geothermal fluid re-injection wells at the site were estimated to have relatively high probability of seismic event occurrence, and these regions experienced high levels of induced seismicity over the microseismic monitoring period. The seismotectonic state estimated for each zone is generally consistent with the forecasted susceptibility to seismicity based on statistics of fracture criticality. A probabilistic seismic susceptibility evaluation framework which integrates seismic observations and hydrogeological modelling is proposed. Fracture criticality is proposed as a novel reservoir depth-independent metric for fracture slip susceptibility. Fracture criticality values estimated at Hellisheiði roughly follow a Gaussian distribution ranging between 0.001 and 2.0 bar/km. Satisfactory probabilistic seismic susceptibility evaluation results were achieved at Hellisheiði. A probabilistic seismic susceptibility evaluation framework which integrates seismic observations and hydrogeological modelling is proposed. Fracture criticality is proposed as a novel reservoir depth-independent metric for fracture slip susceptibility. Fracture criticality values estimated at Hellisheiði roughly follow a Gaussian distribution ranging between 0.001 and 2.0 bar/km. Satisfactory probabilistic seismic susceptibility evaluation results were achieved at Hellisheiði.


Introduction
Induced seismicity is a longstanding concern in geothermal production, where energy-depleted geothermal fluids are re-injected into the subsurface and circulate in the reservoir system.Despite extensive research on the mechanism of induced seismicity (e.g., Rathnaweera et al. 2020), evaluating the susceptibility to fluid injection-induced seismicity remains a significant challenge.
Induced seismicity assessment is a key component to inform seismic hazard evaluation (seismicity that has the potential to cause damage) and seismic risk (the probability of unfavourable consequences of seismic hazards on humans and their built environment) (Schultz et al. 2021).This requires a thorough understanding of conditions under which seismic events are induced by fluid injection.A rigorous analysis of susceptibility to injection-induced seismicity involves careful examination of three different properties: (1) hydrological and geomechanical response of subsurface reservoirs under injection conditions; (2) heterogeneous distribution of pre-existing fractures in the subsurface, which are difficult to directly detect and characterise; and (3) natural variability of fracture attributes and rock properties, which can result in the variability of fluid pressure perturbations required to trigger seismicity.Although either rock failure or fracture slippage in the subsurface is governed by the deterministic Mohr-Coulomb failure/slippage criterion, the intrinsic variability of rock properties and fracture attributes inevitably brings about uncertainty and stochasticity in the seismic generation process.
Hydrogeological modelling targets the first property by representing the deterministic physical process from a mechanics perspective, which forms the basis for physics-based seismic evaluation methods (Zhang et al. 2013;Grigoli et al. 2017).Besides outcrop fracture mapping and geological structure characterisation, observations of microseismicity, which result from the slippage of a subset of underlying pre-existing fractures, are a truthful reflection of the spatial distribution and density of fractures, and contribute to the understanding of the second property (Fisher et al. 2004;Zhao et al. 2019).Probabilistic assessment of induced seismicity proves an effective tool to address the third property by quantifying uncertainties involved in the seismic generation process (Rothert and Shapiro 2007;Walsh III and Zoback 2016;Seithel et al. 2019;Shen et al. 2019).
Having been used for several decades, probabilistic seismic evaluation is emerging as a standard methodology to mitigate fluid injection-induced seismicity in recent years (Grigoli et al. 2017).One of the most prominent probabilistic approaches used is the probabilistic seismic hazard analysis (PSHA), a methodology used to estimate the likelihood of exceeding a certain ground motion intensity over a given time period in the future (Cornell 1968).The PSHA framework typically includes two main components, a seismogenic source model that provides spatiotemporal forecasts of earthquakes, and a ground motion characteristics model that links the earthquake with ground motion.PSHA has been applied to develop earthquake-triggering probability maps in many areas, though its utilisation has constantly been debated and questioned (e.g., Mulargia et al. 2017).
This method assumes long-term stationarity of earthquake rates, therefore, its application to fluid injection-induced seismic risk evaluation may be inherently problematic due to the highly varying injection rates and consequent seismic rates (Langenbruch et al. 2018).Nevertheless, the PSHA has been extended by considering subsurface fluid injection activities to develop regional induced seismic hazard risk estimates (e.g., Gupta and Baker 2019;Broccardo et al. 2020).In addition, a number of risk evaluation methods that use a probabilistic scenario-based event tree approach, similar to the PSHA, were proposed to estimate the overall fault system reliability by combining the likelihood of different scenarios and the failure probability of each scenario (e.g., Mignan et al. 2015;Rahman et al. 2021).
Another group of probabilistic seismic evaluation methods are based on the probabilistic fault stability evaluation by considering the natural variability of fault mechanical properties and stresses (Rothert and Shapiro 2007;Walsh III and Zoback 2016;Seithel et al. 2019;Shen et al. 2019).Besides, the variability inherent to fault mechanical properties and stresses, uncertainties are amplified by two selfenforcing effects during fluid injection processes, i.e., stress field alteration by hydraulic and thermal effects, and degradation of fault attributes by carbonate dissolution (Seithel et al. 2019).Uncertainties considered involve those associated with each Mohr-Coulomb parameter, such as stress, pore pressure, coefficient of friction, and fault orientation.The quantitative susceptibility assessment involves the slip tendency analysis using Monte-Carlo simulations to calculate the cumulative distribution function of the critical pore pressure required to cause slip, based on different combinations of input parameters.In these methods, the fluid pore pressure utilised throughout the investigated region was approximated in different ways, such as being approximated by an analytically derived evolving value assuming diffusion from a point injection source (Rothert and Shapiro 2007), an assumed uniform value based on previous modelling efforts (Rahman et al. 2021), the ambient formation pore pressure value measured from several boreholes (Shen et al. 2019), or a random value drawn from a statistical distribution (uniform distribution, e.g., Walsh III and Zoback 2016;normal distribution, e.g., Seithel et al. 2019).In this respect, recorded seismic events were not accurately correlated with prevailing fluid pressure changes that drive the seismic occurrence, which could lead to errors in estimating the critical fluid pressure.
Recently, efforts have been made to combine probabilistic seismic susceptibility evaluation with hydrogeological modelling, which provides a more accurate spatiotemporally varying fluid pressure distribution.Dempsey and Suckale (2017) carried out a study forecasting induced seismicity due to natural gas extraction using a combination of a physics-based model to account for poroelastic earthquake triggering, and an ensemble approach to represent uncertainties related to the fault strength and stress state, where probability distributions of material and stress parameters are determined from the earthquake data using a Bayesian method.Langenbruch et al. (2018) forecasted the likelihood of damage inducing earthquakes in both space and time using a hybrid physical-statistical model, which involved the evaluation of pore pressure change by a regional hydrogeologic model, and analysis of spatial variability of the seismogenic state.Hennings et al. (2021) integrated the pore pressure change derived from hydrogeological modelling, earthquake catalogues, and probabilistic fault slip potential analysis to investigate relationships between earthquake sequences and the spatiotemporal association with oilfield wastewater disposal activities.In addition, integration of geomechanical analysis and microseismic observations has been used for various evaluation purposes, such as to estimate the fracture pressure to ensure CO 2 storage reservoir integrity (Goertz-Allmann et al. 2014), and to map seismically active regional fault-containing corridors (Weir et al. 2022).However, there are only a few induced seismicity evaluation methods that incorporate seismic observations, hydrogeological modelling and probabilistic seismicity occurrence analysis.Considering that the fracture slippage criterion is satisfied at the source fracture for each seismic event, the integrated analysis of seismic observations and hydrogeological modelling allows for quantifying the variability of perturbations required to trigger seismicity, thus providing a probabilistic perspective in evaluating the susceptibility to injection-induced seismicity.
This paper presents work carried out on the development and application of a probabilistic seismic susceptibility evaluation framework by combining seismic observations and hydrogeological modelling of fluid injection operations in faulted reservoirs.Fracture criticality, defined as the gradient of critical fluid pressure change to trigger seismicity (Δp c /h), is proposed as a novel reservoir depthindependent metric of fracture slip susceptibility and used in the probabilistic evaluation.This evaluation framework is demonstrated using induced seismicity data recorded during a 2-year fluid re-injection period (15 April 2020-15 April 2022) at the Hellisheiði geothermal field, Iceland.During the probabilistic evaluation, induced seismicity data in a half-year period were first analysed and spatially correlated with major fault structures.Hydrogeological simulations were then conducted to represent the hydrological response to the fluid injection over the re-injection period.The integrated interpretation of fracture criticality was performed for each fault zone, and was used to provide probabilistic injection-induced seismic hazard maps for the next half-year period at the site.

A Seismic Susceptibility Evaluation Method Based on Statistics of Fracture Criticality
Slippage of pre-existing fractures has been recognised as the dominant mechanism for fluid injection-induced seismicity (e.g., Raleigh 1972;Dieterich 1978;McGarr 2014;Elsworth et al. 2016).In this respect, the susceptibility to fluid injection-induced seismicity can be evaluated in terms of the fracture slip susceptibility.For any selected measure of fracture slip susceptibility, the slippage of each individual pre-existing fracture is characterised by a critical value, according to the Mohr-Coulomb fracture slippage criterion.When this measure of fracture slip susceptibility is independent of reservoir depth, the statistical distribution of its critical values for a set of slipped fractures represents only uncertainties inherent in the fracture slippage process, which could be used for probabilistic forecasting of future seismicity.
Based upon this understanding, a probabilistic seismic susceptibility evaluation method that integrates recorded seismicity interpretation and hydrogeological modelling is developed to estimate the seismicity occurrence probability, given that (1) statistics of fracture slip susceptibility, (2) the density of underlying fractures, and (3) prevailing stress and fluid pressure conditions, are known.As a prerequisite to this method, a suitable metric of fracture slip susceptibility that is independent of the reservoir depth needs to be selected.Recorded seismicity can be used as a gauge for the critical fracture slip susceptibility for the probabilistic evaluation (Grasso et al. 1992) and a measure of the regional fracture density (Cao et al. 2020a, b).Hydrodynamic modelling of fluid injection into the reservoir can be used to obtain stress and fluid pressure distribution at the time of forecasting.
Section 2.1 summarises conventional metrics of fracture slip susceptibility based on the Mohr-Coulomb fracture slippage criterion, followed by two novel reservoir depthindependent metrics suitable for probabilistic seismic susceptibility evaluation in Sect.2.2.In Sect.2.3, a probabilistic approach to evaluate induced seismicity based on statistics of fracture slip susceptibility is discussed.The procedure which combines seismic observations and hydrogeological modelling in evaluating probabilistic seismic susceptibility is presented.

Conventional Metrics of Fracture Slip Susceptibility
Shear reactivation of an individual fracture can be described by the Mohr-Coulomb fracture slippage criterion: (1) where τ is the shear stress along the fracture, σ and p are the normal stress and fluid pressure applied on the fracture plane, and f and c are the friction coefficient and cohesion of the fracture, respectively.For pre-existing fractures, c ≈ 0 , and Eq.(1) can then be simplified to: The normal and shear stresses are given by: where σ 1 and σ 3 are the maximum and minimum principal stresses, respectively, and β is the angle between the outward unit normal vector of the fracture and the direction of σ 1 .It is noted that this criterion is different from the Mohr-Coulomb failure criterion for intact rocks, where the most vulnerable orientation to rupture makes an angle of (45°-φ/2) with the maximum principal stress orientation, φ being the friction angle.
The Mohr-Coulomb fracture slippage criterion Eq. (2) forms the basis for a number of indices that have been widely used to evaluate the critical state of pre-existing fractures and the associated susceptibility to seismicity under fluid injection.These indices include: (1) fluid pressure p (e.g., Rothert and Shapiro 2007;Chiaramonte et al. 2008;Shapiro 2015;Kettlety et al. 2021;Shen et al. 2021), (2) fluid pressure change Δp (e.g., Davis and Pennington 1989;Rajendran 1995;Hainzl et al. 2012;Shen et al. 2021), (3) Coulomb failure stress change ΔCFS (defined as the static stress change associated with displacement along a fracture plane) (e.g., Cao et al., 2021Cao et al., , 2020b;;King et al., 1994;Reasenberg and Simpson, 1992;Scholz, 2002), and (4) Slip tendency T S (defined as the ratio of shear stress to normal stress on a fracture plane) (e.g., Morris et al. 1996;Blöcher et al. 2018;Shen et al. 2021).Table 1 lists some conventional metrics of fracture slippage susceptibility based on the Mohr-Coulomb fracture slippage criterion, and their expressions and applicability.Amongst these indices, the Coulomb failure stress change ΔCFS and slip tendency T S reflect the fracture slip susceptibility subjected not only to the fluid pressure change, but also perturbing stresses on the receiving fracture surface.Both indices are suitable to evaluate susceptibility to seismicity induced by a broad range of triggering mechanisms, such as rising fluid pressure (Talwani and Acree 1985;Elsworth et al. 2016), poroelastic stressing (Segall 1989;Segall and Lu 2015), thermoelastic stressing (Segall and Fitzgerald 1998), and fault slippage-induced stress transfer (Guglielmi et al. 2015;Cao et al. 2021).In contrast, fluid pressure p and fluid pressure change Δp can be used to evaluate potential  for seismicity caused predominantly by the fluid overpressure effect, when normal and shear stress changes on preexisting fractures are negligible.In this case, the Coulomb failure stress change ΔCFS is reduced to be proportional to the fluid pressure change Δp: ΔCFS = fΔp.Critical values of these indices required to trigger seismicity can be derived directly from the Mohr-Coulomb fracture slippage criterion Eq. (2).As shown in Table 1, p c , Δp c , ΔCFS c are expressed as a function of prevailing normal and shear stresses σ and τ and the friction coefficient f, whilst T Sc equals to the friction coefficient f.As such, statistics of critical values for the first three indices account for uncertainties related to both tectonic stresses and friction coefficient (King and Cocco 2001;Scholz 2002), whilst those for the slip tendency, T Sc only reflect uncertainties relevant to the friction coefficient.
Figure 1 graphically illustrates critical values required to trigger seismicity for the four indices.In situ stresses on confined rock masses are represented by blue Mohr's circles, whilst the in situ fluid pressure p shifts the stress states to orange Mohr's circles approaching the Mohr-Coulomb fracture envelope.Considering a pre-existing fracture plane whose outward unit normal vector makes an angle of β with the maximum principal stress orientation, its stress states and effective stress states are, respectively, represented by points A and B on Mohr's circles.p c , Δp c , ΔCFS c and T Sc are denoted by the distance from point A horizontally to the Mohr-Coulomb envelope (AC), distance from point B horizontally to the Mohr-Coulomb envelope (BC), distance from point B vertically to the Mohr-Coulomb envelope (BD), and the slope of OB, respectively.Herein, the critical fluid pressure p c (= p + Δp c ) required to trigger seismicity is defined as rock criticality (Rothert and Shapiro 2007;Shapiro 2015).It is noteworthy that p c , Δp c and ΔCFS c generally increase with the reservoir depth, as such they may vary markedly across the depth range of field monitored seismicity (Fang et al. 2016).

Fracture Criticality as a Reservoir Depth-Independent Metric of Fracture Slip Susceptibility
To account for uncertainties intrinsic to the fracture slip process, the fracture slip susceptibility index used for probabilistic evaluation of induced seismicity needs to meet the following requirements: (1) critical value to trigger fracture slippage is independent of reservoir depth, and (2) statistics of the critical value reflects uncertainties in both tectonic stress gradients and fracture frictional properties.As such, none of the four indices discussed above meets both requirements.The fluid pressure p, fluid pressure change Δp and the Coulomb failure stress change ΔCFS, being generally depth dependent, fail to meet the first requirement, whereas the slip tendency T S fails to meet the second requirement, as statistics of its critical values would only involve uncertainties in fracture frictional properties, but not in tectonic stresses.To this end, two novel reservoir depth-independent metrics of fracture slip susceptibility, termed the gradient of Coulomb failure stress change ΔCFS/h and the gradient of fluid pressure change Δp/h, are proposed in this work.When normal and shear stress changes on pre-existing fractures are negligible, the two indices are related by a linear relationship: ΔCFS/h = fΔp/h.In many applications, in situ stresses and fluid pressure values (σ, τ and p) may be approximated as linear functions of reservoir depth h, in the form σ = a 1 h, τ = a 2 h, and p = a 3 h (where a 1 , a 2 and a 3 are constants).By substituting these relationships to expressions for the critical Coulomb failure stress change ΔCFS c and critical fluid pressure change Δp c (see Table 1), critical values for the two indices can be expressed as: (5) It can be seen that these two indices meet both the requirements and therefore are able to characterise the perturbations required to trigger seismicity (or the critical state of pre-existing fractures).The lower is ΔCFS c /h or Δp c /h, the more critical the pre-existing fracture.
ΔCFS c /h or Δp c /h for induced seismicity recorded (or pre-existing fractures that undergo slippage) is also easy to obtain by integrating microseismic observations and hydrogeological simulation as discussed later.The gradient of fluid pressure change Δp/h is used as the metric of fracture slip susceptibility in probabilistic seismic susceptibility evaluation in this work.For convenience, the gradient of critical fluid pressure change Δp c /h is defined as fracture criticality C.

Seismic Susceptibility Evaluation Based on Statistics of Fracture Criticality
Owing to the inherent heterogeneity of tectonic stresses, fracture attributes (such as fracture orientation) and rock properties (such as friction coefficient), the criticality (C = Δp c /h) of a random set of fractures subjected to fluid injection may vary in a range (C min ~ C max ), even within the same host rock.The probability for an individual preexisting fracture to satisfy Mohr-Coulomb fracture slippage criterion equals to the probability of the prevailing gradient of fluid pressure change (Δp/h) being greater than the fracture criticality.Considering that the fracture criticality is statistically distributed, the probability can be expressed as (Shapiro 2015): where f(C) is the probability density function (PDF) of the fracture criticality C(r) at any spatial location r.Given the observed range for the fracture criticality (C min ~ C max ), the seismic occurrence probability would be 0 when Δp/h < C min , and 1 when Δp/h > C max .
In the case of the simplest possible PDF for the fracture criticality, i.e., the uniform PDF where , the seismic occurrence probability can be simplified as a linear function of Δp/h: Considering the evaluation of seismic susceptibility in one gridblock, the probability of seismic event occurrence within this particular gridblock region can be expressed as a function of the probability for an individual fracture to satisfy the Mohr-Coulomb fracture slippage criterion P f , and the regional fracture density d within the gridblock region: When the regional fracture density d ≤ 1 per gridblock, the probability of seismic event occurrence is directly calculated by the product of P f and d.When the regional fracture density d > 1 per gridblock, the probability of seismic event occurrence is calculated based on the probability that no seismicity occurs (1 − P f ) d .
The injection-induced seismic susceptibility evaluation described above can be performed for injection field sites containing various geological structures or stratigraphical formations.Key assumptions of the seismic evaluation method include: (1) a random set of pre-existing fractures as hypocentres of induced seismicity is considered to be statistically homogeneously distributed within each geological structure or stratigraphical formation; (2) pre-existing fractures within different geological structures or stratigraphical formations are characterised by different fracture attributes (such as fracture density and fracture criticality); (3) fractures are orientated parallel to the regional fault strike, and the fracture density is proportional to the density of field monitored seismicity (Cao et al. 2020a, b); and (4) pre-existing fractures do not mutually interact.
The evaluation starts with the zonation of geological structures in the region of investigation.The areal extent of geological structures (e.g., fault zones) is identified based on available data, including structure models for fault systems, outcrop fracture mapping, and historical seismic clustering observations.From the start of fluid injection operations, a sequence of fluid injection parameters (such as injection rate and wellhead pressure) and concurrent seismic catalogue are recorded and divided into sequential periods.The evaluation involves integrated seismic and hydrogeological analysis for each period, and probabilistic forecasting of seismicity for the next period.This evaluation process, besides the comparison between recorded seismicity and probability forecasts, is iterated for each period until the end of fluid injection operations.
The procedure for probabilistic forecasting through integration of microseismic observations and hydrogeological simulation comprises the following steps: Step 1, analysis of field monitored induced seismicity within the period considered is performed.3D hypocentre coordinates of recorded seismicity are used to filter ones falling within the region corresponding to the reservoir model domain.Each recorded seismic event within a specified distance to known geological structures is assigned to the nearest one.The fracture density in each geological structure is inferred from the density of induced seismicity falling within this zone.( 9) Step 2, hydrogeological modelling of reservoir response to fluid injection is carried out.To ensure the reservoir behaviour to be realistically represented, the reservoir model constructed would incorporate known geological structures and be calibrated through historically matching field injection pressure time series over the injection period according to field injection rates.The fluid pressure change Δp at the hypocentre of each recorded seismic event at the time of seismic occurrence is considered as the critical fluid pressure change Δp c to trigger the event, and is extracted from the calibrated reservoir model.The fracture criticality Δp c /h is then calculated for each induced seismic event.
Step 3, estimation of the susceptibility to seismicity occurrence based on the statistics of fracture criticality is achieved.The histograms of the fracture critical-ity, representing its probability density, can be obtained for each individual geological structure, or collectively by aggregating all possible values.Based on modelled fluid pressure change Δp throughout the model domain for the next period, the likelihood that the Mohr-Coulomb fracture slippage criterion is satisfied at any spatial location for the next period can be calculated according to Eq. ( 7).Combining the fracture density obtained in Step 1, the likelihood of seismic event occurrence within each gridblock throughout the model domain for the next period is estimated according to Eq. ( 9).A flowchart of the procedure for the seismicity occurrence probability forecasting is presented in Fig. 2. Step 3: Seismic susceptibility evaluation

Geothermal Fluid Re-injection at Hellisheiði Geothermal Field
Located in the Southwest Iceland, the Hellisheiði geothermal field is one of the largest geothermal fields in the country, producing geothermal fluid at 240-320 °C from a subsurface reservoir at depths ranging from 1.5 to 3.3 km.The Hellisheiði geothermal heat and power plant, operated by Orkuveita Reykjavíkur (Reykjavik Energy), started production in 2006, and the current installed capacity is 303 MWe for power generation and 200 MWth for space heating.
The general stratigraphy at Hellisheiði comprises of alternating successions of hyaloclastite formations formed during glacial periods and lava sequences from interglacial periods.The region has a complicated tectonic setting, consisting of NNE-SSW trending extensional fault structures with sub-vertical dips and highly fractured fissure swarms (Snaebjörnsdóttir et al. 2018;Ratouis et al. 2019) (Fig. 3).These faults represent the western margin of the Hengill fissure zone (Batir 2011).Faults in the Hellisheiði region are normal faults, whilst faults in the South Iceland Seismic Zone are mostly strike slip faults with N-S orientation (Gunnarsson 2013).The maximum horizontal principal stress follows a NNE orientation, aligned with the strike of major faults (Batir 2011).
Operation permits require all geothermal brine after production to be re-injected back to the geothermal reservoir.Re-injection of the geothermal brine and condensate takes places in two distinct sites.The main re-injection site, the Húsmúli area, is located on the north-western edge of the field, and annually re-injects around 12 Mt of geothermal spent fluid through five boreholes (HN09, HN12, HN14, HN16 and HN17).The re-injection wells at Húsmúli target major fault structures: the fault F1 cuts across well trajectories of HN14 and HN17, the fault F3 cuts across those of HN09, HN14, HN16 and HN17, and the fault F4 cuts across those of HN09, HN12, HN16 and HN17 (Fig. 3).Wells are cased off down to ~ 800 m depth to obstruct the loss of circulation fluid (Snaebjörnsdóttir et al. 2018).Fluid re-injection into wells HN09, HN12, HN14 and HN17 started on 1 September 2011, followed by the well HN16 on 23 September (Juncu et al. 2020).Since 2014, spent geothermal fluid reinjected into well HN16 has been dissolved with produced CO 2 to reduce carbon emissions from the Hellisheiði power plant.
Continuous microseismic monitoring has been carried out by a regularly extended regional seismic network over the Hengill area since 1993.The commissioning of the Húsmúli re-injection site has been accompanied by swarms of induced seismicity.Whilst a number of induced seismicity emerged during the drilling and pumping tests, the intensity of induced seismicity became unprecedently high shortly after the operation of the site and faded out over time (Gunnarsson 2013).Historical microseismic observations have shown that induced seismicity tends to occur at faults cutting the re-injection wells at the site (Hjörleifsdóttir et al. 2019).

Seismic Susceptibility Evaluation Based on Statistics of Fracture Criticality at Hellisheiði Geothermal Field
Fluid injection operational data and concurrent seismic catalogue during a 2-year period between 15 April 2020 and 15 April 2022 at the Hellisheiði geothermal field were used to demonstrate the induced seismicity susceptibility evaluation method based on statistics of fracture criticality.
The seismic catalogue was kindly provided to the authors by the COSEISMIQ project.This seismic catalogue was divided into four sequential half-year periods, such that integrated microseismic interpretation and hydrogeological modelling over early periods could be used to provide seismic occurrence forecasts over subsequent ones.Sections 4.1-4.3showcase the field microseismic data analysis, hydrogeological modelling of concurrent geothermal fluid re-injection, and probabilistic seismic evaluation for the first half-year period (between 15 April and 15 October 2020), which provides seismic occurrence forecasts over the second half-year period (between 16 October 2020 and 15 April 2021).Section 4.4 presents the half-year seismicity occurrence probability forecasts for four half-year periods and comparison against recorded seismicity when seismic records are available.

Geothermal Fluid Re-injection-Induced Seismicity
Figure 4 presents field recorded operational time series and monitored seismicity for geothermal fluid injection into the five re-injection wells during the first half-year period at Húsmúli area.With wellhead pressures being controlled at ~ 8 bar, the five re-injection wells can be considered as constant pressure boundaries.From field monitored injection operational time series, injection rates and wellhead pressures at the five re-injection wells are relatively stable (top four panels in Fig. 4).As a result, the field well injectivity (defined as the fluid injection rate over the wellhead pressure) remains fairly stable over this period.The fluid circulation through the geothermal system can be considered to have reached quasi-steady-state flow conditions, under which the spatial distribution of the fluid pressure in the reservoir remains almost constant over time.
Induced seismic activities were persistently occurring at Húsmúli over the period of investigation, totalling at 975 recorded seismic events (the bottom panel in Fig. 4).Whilst daily event counts mostly remain below 20, seismic activities with daily event counts above 40 were recorded in late April, early August and early September 2020.Daily maximum seismic magnitude M L ranges between 0 and 3.No conspicuous correlations have been identified between daily seismic event counts, daily maximum seismic magnitude and fluid injection operations.
Figure 5 presents the spatial distribution of induced seismicity over the period of investigation at the Húsmúli area.The depth of induced seismic events shows significant variation, ranging from close to the ground surface down to over 6 km depth.A large proportion of the seismic events are located in the vicinity of major fault structures.Whilst induced seismicity went through a quiescence period during 16 May-15 July (Fig. 5b, c), they were relatively active over the remaining period (Fig. 5a, d-f).The largest seismicity swarm generally located in the vicinity of major faults, shifting from the F1 fault zone during 15 April-15 May (Fig. 5a), to the F2 fault zone during 16 July-15 September (Fig. 5d,  e), and back to the F1 fault zone during 16 September-15 October (Fig. 5f).This marked volatility in spatial and temporal seismic characteristics were also reflected in historical seismic observations between September 2011 and May 2012 at the same field site (Juncu et al. 2020).Nevertheless, it is difficult to identify the respective contribution in triggering these seismic events from re-injection into each well due to the proximity between well trajectories of the five injection wells.
Out of all the events recorded during 15 April and 15 October 2020 at Húsmúli, a total of 718 seismic events were identified to locate within the reservoir model domain as discussed later in the following Sections.Each seismic event within 200 m distance to the major faults was assigned to the closest fault.The remaining seismic events (those located at off-fault zones) were categorised as a separate group.Whilst 623 seismic events (86.8%) were associated with faults, 95 events (13.2%) located at off-fault zones.

Hydrogeological Modelling of Geothermal Fluid Re-injection at Hellisheiði
A reservoir model of a part of the Húsmúli area measuring 6 km × 5 km × 3 km (northing length × easting length × height) was constructed in the industry standard reservoir simulator ECLIPSE 300 to simulate the reservoir response of the five major faults and surrounding areas to fluid re-injection from the five re-injection wells over the period 15 April-15 October 2020.The gridblock dimension of the model is 100 m in all the three directions.Five fault structures (F1-F5), which represent planes of crustal weakness and fracturing resulting from intrusions and fissure eruptions, were incorporated in the reservoir model (Fig. 6).Digitised from the geological model, these faults were extended to the surface to fit the surface fault map, and down to 3 km depth to allow fluid flow at varying depths.The fault F2 was identified as a possible fault by Snaebjörnsdóttir et al. (2018) and further confirmed by Cai et al. (2021) from the interpretation of seismic monitoring results in the area.F2 was implemented in the model as vertical planes extending from surface lineaments.
Both the reservoir and basement formations in the model domain were assumed to have uniform reservoir properties (porosity and permeability) for simplicity.According to the analysis carried out on 80 samples of basaltic hyaloclastic tuffs (Frolova et al. 2004;Snaebjörnsdóttir et al. 2014), the hyaloclastic formation as a reservoir manifests broad variations in both porosity and permeability, varying in the range of 0.14-0.57and 0.001-6400 mD, respectively.Snaebjörnsdóttir et al. ( 2018) used a porosity of 0.10 and a permeability of 30 mD for the hyaloclastic formation in their hydrogeological modelling study.In this work, the porosity value of the reservoir formation was set to 0.10, while the permeability values for the reservoir formation and the faults were determined from history matching the bottomhole pressure time series over the period of investigation.The combination of a permeability of 60 mD for the reservoir formation, and a much higher permeability of 60 Darcy for the faults were found to provide a reasonable match for the five reinjection wells (Fig. 7).The basalt basement was assigned a porosity of 0.08 and a permeability of 3.75 mD, following Stefánsson et al. (1997) and Snaebjörnsdóttir et al. (2018).A pore volume multiplier of 10,000 was applied to the lateral boundaries of the reservoir model to account for the large volume of pore space surrounding the model domain.
The natural groundwater table at the Húsmúli area is ~ 200 m above the sea level, which results in an overpressure of 20 bar in the geothermal reservoir (Gunnarsson 2013).The geothermal reservoir model was first initialised to hydrostatic conditions with this designated overpressure (Fig. 6).According to field injection rate and pressure time series (Fig. 7), fluid circulation in the geothermal reservoir at the Húsmúli area had already been established before the modelling started, and fluid flow has maintained a quasisteady-state condition throughout the modelling period.Therefore, the hydrological effect from the previous decade of fluid injection since 2011 is considered to have limited influence on the model results.Nevertheless, a pre-modelling step, accounting for fluid injection over a month's period before 15 April 2020, was then simulated so as to establish fluid circulation in the geothermal system.
The history matching efforts attempted to match bottomhole pressures for all five injection wells at the same time.Figure 7 presents the history matching results for the five re-injection wells at the Húsmúli area.The bottomhole pressures were calculated from the field monitored wellhead pressures, with consideration of the overpressure in addition to hydrostatic pressure at the reservoir depth.The modelled and field injection bottomhole pressures achieve a fairly close match over the period of investigation, indicating that the fluid pressure field obtained from this reservoir model could represent the critical fluid pressure to induce seismicity at the Húsmúli area.It is noted that pressure time series for well HN16 are characterised by more fluctuations than others, likely due to the concurrent CO 2 injection at this well, which was not considered in the hydrogeological modelling.Sources of deviations in history matching for other wells include heterogeneous distribution of porosity and permeability, and fluid injection and thermal contraction-induced permeability change, which were not explicitly modelled in this work.Figure 9 presents probabilistic seismic susceptibility estimation results for the second half-year re-injection period at 1 km depth at the Húsmúli area.As shown in Fig. 9a, the modelled re-injection-induced fluid pressure increase is around 2 bar surrounding the five re-injection wells, and it decreases away from the wells.Largely dictated by the major faults, the contours of fluid pressure change are predominantly orientated in alignment with the fault strike.There are also local extensions of the contour boundary around the faults, resulting from larger elevated fluid pressure changes due to the higher permeability of faults.Owing to quasisteady-state fluid flow conditions reached at the Húsmúli area, the spatial distribution of fluid injection-induced pressure change is maintained almost constant over the period of investigation.The probability of satisfying the fracture slippage criterion (Fig. 9b) is influenced by both the fluid pressure increase driven by fluid injection operations (Fig. 9a), and the intrinsic variability of fracture criticality related to the natural variability of fracture attributes and rock properties (Fig. 8).This probability reaches close to 100% within a large proportion of the model domain, including F3 and F4 fault zones with a Δp/h greater than 1.0 bar/km, F1, F2 and F5 fault zones with a Δp/h greater than 0.5 bar/km, and offfault zones with a Δp/h greater than 0.3 bar/km (Fig. 9a).Influenced by both the likelihood of satisfying the fracture slippage criterion (Fig. 9b) and regional fracture density related to neighbouring major faults (Fig. 9c), relatively high probability of seismic event occurrence was estimated for fault zones around the five re-injection wells (Fig. 9d).The fracture density in off-fault zones was estimated as fairly low owing to the much less active induced seismicity in these zones during the first half-year period, therefore, the forecasted seismic susceptibility for these zones during the second half-year period was close to 0.

Comparison of Half-Year Seismicity Occurrence Probability Forecasts Against Recorded Seismicity
Figures 10 and 11 present the forecasted likelihood of seismicity occurrence within each gridblock for four subsequent half-year injection periods at 1 km reservoir depth and within the 3D model domain, respectively.Field monitored seismic events during the first 3 half-year periods are plotted on both graphs for comparison.According to the forecasts, seismically active regions extend from fault zones F1 and F2 during October 2020-April 2021, to fault zone F3 during April-October 2021, and gradually shrink to fault zones F3 and F5 during April-October 2022.The maximum forecasted seismicity occurrence likelihood is estimated to be reached during the period April-October 2021.Assuming that the fluid injection operational conditions remain unchanged, the forecasts suggest a seismic quiescence period during April-October 2022, owing to a much reduced seismic event count recorded (representing the intensity of potentially vulnerable underlying fractures) during the previous half-year period.A plausible explanation accounting for this seismic quiescence period is the Kaiser effect, meaning that seismic events would not occur again unless the prior fluid pressure in the reservoir is exceeded.
Overall, the areas estimated to have high seismic susceptibility are fairly consistent with the seismically active areas during the periods of investigation.The majority of induced seismicity is observed to cluster within fault zones neighbouring the re-injection wells, whilst much less seismic events are sparsely distributed in off-fault zones and fade off away from the re-injection wells.During October 2020-April 2021, most induced seismicity is observed to fall within regions with an elevated fluid pressure in excess of 0.2 bar, in comparison with Fig. 9a.The period October 2020-April 2021 observes relatively intensive seismicity located within off-fault zones, resulting in a distinct elevation in the forecasted seismicity occurrence likelihood within these zones during April-October 2021.It is noted that the recorded seismicity is predominantly located at the northeast end of the forecasted seismically active regions, with few occurring at the southwest end.This deviation indicates that the fault zones assumed to be homogenous in terms of hydraulic conductivities are most likely heterogeneous, forming fluid pathways from the re-injection wells towards the northeast direction but not the southwest direction.
Figure 11 presents plots of regions with elevated seismicity occurrence probability (> 10%) for each half-year forecasting period.The forecasted seismically active regions extend from the ground surface to below 3 km depth, but are characterised by the inverted cone shape.This indicates that regions at shallower depths are more susceptible to seismicity under the same fluid pressure perturbations (see Fig. 1).The fault zones F1 and F2 across the full depth range have remarkably large forecasted seismicity occurrence probability during the first three forecasting periods, whilst seismic activities are forecasted to take place mostly within shallow regions of fault zones F3 and F5 during the fourth forecasting period.Amongst the four forecasting periods examined, the period April-October 2021 has the largest forecasted seismically active off-fault regions.During the first three injection periods of forecasts when the seismic catalogue is available, the majority of recorded seismicity fall within the forecasted seismically active regions, with only a small amount occurring outside these regions.

Implications of Fracture Criticality for the Seismotectonic State of Injection Sites
The fracture criticality concept proposed in this work complements conventional metrics in quantifying fracture slip susceptibility of injection sites.The correlation of fracture criticality with other metrics (Sect.2.2) allows for quantitative comparison between evaluation results using different metrics at certain conditions.Comparison can be made in term of critical values to trigger seismicity.When the overpressure is the dominant seismicity triggering mechanism, fracture criticality is related to the critical Coulomb failure stress change by fh .Considering a friction coefficient of 0.6 for reservoir rocks at Hellisheiði (Arnadóttir et al. 2003;Juncu et al. 2020), ΔCFS c at Hellisheiði is mostly in the order of a few hundredth to a tenth MPa, comparable to critical values to trigger seismicity reported at other injection sites (e.g., Lim et al. 2020;Cao et al. 2021;Kettlety and Verdon 2021).Comparison can also be made in terms of the magnitude distribution of critical values to trigger seismicity.Fracture criticality obtained for the Hellisheiði geothermal site follows the Gaussian distribution within each region, ranging from approximately 0.001 bar/km to around 2 bar/ km (Fig. 8).In contrast, rock criticality reported for both the Soultz-sous-Forêts and Fenton Hill Hot Dry Rock sites roughly follows a uniform or a truncated Gaussian probability density function (Rothert and Shapiro 2007).Rock criticality values distribute in a broad range between 0.001 MPa and approximately 3 MPa for the first site and between 0.001 MPa and approximately 1 MPa for the second, with the upper bounds being several orders of magnitude larger than the lower bounds.In the Fort Worth Basin, USA, the critical fluid pressure change to reactivate earthquake sequences was reported to be characterised by a positively skewed distribution ranging between 0 and 1 MPa with a mean of around 0.05 MPa (Hennings et al. 2021).This different shape of magnitude distribution of critical values to trigger seismicity may be explained by the seismic event depth considered, and the difference and heterogeneity in fracture attributes and stress conditions at different field sites.
The magnitude distribution of fracture criticality reflects both intrinsic fracture friction heterogeneity and in situ stress variability (Eq.6).In cases where the in situ stress can be considered as homogenous, the statistics of fracture criticality represents the variation of friction coefficient of underlying fractures (Zhang and Ma 2020).This indicates that the statistics of fracture criticality may be utilised to constrain the estimation of the friction coefficient of underlying fractures, given sufficient knowledge on in situ stress conditions and injection operational parameters.In practice, it is noted that similar to the range of rock criticality (Shapiro 2015), the range of fracture criticality obtained may be influenced by the extent of fluid pressure increase at locations of induced seismicity.The upper limit of the fracture criticality obtained is probably constrained by the maximum prevalent fluid pressure increase, therefore, one cannot exclude that the actual fracture criticality can be even higher.In addition, regions with high fluid pressure (such as F3 and F4 fault zones) tend to have large fracture criticality, as opposed to small fracture criticality at off-fault zones with much less fluid pressure perturbations (Fig. 8).
Evaluation of fracture criticality combined with recorded seismicity also allows for quantifying the seismotectonic state of injection sites, which is a measure of the level of seismic activity at injection sites when subjected to certain stress or pressure perturbations (Dinske and Shapiro 2013).In this work, the seismotectonic state Ω is defined as the intensity of induced seismicity N u (the seismic event count per unit area) divided by the fracture criticality (the critical value to trigger seismicity), i.e., Ω = N u Δp c ∕h .In comparison, the seismotectonic state is conventionally represented by the seismogenic index Σ, given by Σ = a + log 10 ( ) , where a = log 10 (N M L ≥0 ) is a seismotectonic constant in the Gutenberg-Richter scaling law, N f is the concentration of fractures, S is the poroelastic uniaxial storage coefficient, and p c,max is the maximum value of rock criticality (Shapiro et al. 2010).Herein, the quantity F t = p c,max /N f is referred to as the tectonic potential, which depends on tectonic activities of the injection region only.The seismic event count N is related to the concentration of fractures N f by N = , where Q c is the cumulative injected fluid volume.By these definitions, it is straightforward that the seismotectonic state Ω at a given reservoir depth h is in general positively correlated with the seismogenic index Σ.
To evaluate the seismotectonic state of different regions at Húsmúli, recorded seismic counts against the average fracture criticality at different regions over time are plotted in Fig. 12a.The seismotectonic state of each region is represented by the slope of the connection line between the corresponding scatter point and the original point, i.e., the steeper the slope, the more susceptible seismogenic conditions the region has.The evolution of the seismotectonic state of different regions over the four half-year injection periods of investigation is presented in Fig. 12b.It can be seen that the seismotectonic states vary to a large extent depending on the region and the time.During the four periods, the seismotectonic state of F1 and F2 fault zones is the most critical, followed by that of F3 and F5 fault zones, whilst that of the F4 fault zone is relatively lower.In comparison, seismotectonic state of off-fault zones is the least critical.Whilst F1, F2 and F3 fault zones evolved seismotectonically from an active state during the first two half-year periods to a much suppressed state during the last two half-year periods, the seismotectonic state of F4 and F5 fault zones and off-fault zones maintained fairly low levels throughout.Evolutionary trends of the seismotectonic state Ω of different regions (Fig. 12b) are in general consistent with those of the seismogenic index Σ (Fig. 12c), with deviations likely caused by the use of different fracture slip susceptibility metrics (fracture criticality Δp c /h against rock criticality p c ).
It is noteworthy that, although the fracture criticality of off-fault zones may be lower than that of faulted zones (fractures are more susceptible to fluid pressure changes) (Fig. 8), the fairly low density of seismicity observed at off-fault zones (Figs.9c and 12a) leads to an overall less critical seismotectonic state of off-fault zones (Fig. 12b) and consequently low estimated susceptibility to seismicity (Fig. 9d).The seismotectonic state estimated for each region during each half-year modelling period (Fig. 12b) is generally consistent with the seismic susceptibility forecasted based on the statistics of fracture criticality (Fig. 10).

Limitations and Applicability of the Probabilistic Seismic Susceptibility Evaluation Method
The seismic susceptibility evaluation method presented here is also a physics-based method, as it has explicitly considered the hydrological behaviour of a reservoir in response to fluid injection.In this regard, accurately representing the spatial and temporal evolution of fluid pressure change within the reservoir is crucial to achieving satisfactory forecasting performance.The work presented here implemented five fault structures in the reservoir model developed based on regional geology and historical seismic observations.Uniformly homogeneous hydrological properties were distributed within each region so that the entire fault structures were considered as preferential fluid flow paths.This may have led to less satisfactory forecasting results at some localised regions, such as the southwest end of the forecasted seismically active regions within which few recorded seismicity falls (Fig. 10).Nevertheless, anisotropy in permeability distribution at the reservoir level has been introduced to reflect the tectonic controls on fluid flow paths by permeabilities of the five major faults in the geological model.In comparison, previous modelling studies have represented fluid flow paths at Húsmúli using a 1D flow channel representation (Kristjánsson et al. 2016) 2021) inferred three preferential fluid flow paths corresponding to three major NNE-SSW trending faults in the geological model.A constant fracture permeability of 87.5 mD was assigned for the fluid flow paths (assuming an average thickness of 120 m of gridblocks), whilst the remaining model domain had a fracture permeability of 0.5 mD.The calibrated model implementing these fluid flow paths could match results of an independent tracer test with acceptable accuracy.In this regard, the forecasting performance of the seismic evaluation method may be improved by considering these finer scale preferential fluid flow paths in hydrogeological modelling.
The seismic susceptibility evaluation method proposed was used to evaluate and forecast seismicity driven by the direct increase in fluid pressure, which is the most common seismicity triggering mechanism.Here, the gradient of critical fluid pressure change Δp c /h is utilised as the fracture slip susceptibility metric.However, it is acknowledged by the authors that in many subsurface fluid injection applications, particularly in geothermal systems where cold fluids are injected to reservoirs, a number of seismicity triggering mechanisms are at play, including the elevated fluid pressure, poroelastic stress, thermal effects, interactions between seismicity, and the Kaiser effect (Izadi and Elsworth 2015;Ghassemi and Tao 2016;Brown and Ge 2018;Garcia-Aristizabal 2018;Im and Avouac 2021;Cao et al. 2022).The overpressure-driven seismicity has a larger areal extent of influence and occurs rapidly after injection, whilst the temperature-driven seismicity is more constrained around injection wells and may be dominant in the long term (Izadi and Elsworth 2015;Ghassemi and Tao 2016).The direct pore pressure change enhances the potential for seismicity in overpressurised regions, whilst the poroelastic stressing reduces the potential for seismicity in overpressurised regions but enhances that surrounding these regions (Cao et al. 2022).The thermal stress is dominant in the neighbourhood of injection wells, enhancing the potential for seismicity in regions perpendicular to the maximum horizontal stress direction, and inhibiting that in regions along the maximum horizontal stress direction (Cao et al. 2022).The relative contribution of the hydrological and thermal effects depends on the distance from seismicity to injection wells and the injection duration.The seismicity recorded during the modelling period may be due to thermal effects from the previous decade of fluid injection, in addition to the direct pore pressure effect.When considering thermal effects for more accurate probabilistic seismic evaluation, the gradient of critical Coulomb failure stress change ΔCFS c /h can be used as the metric of fracture slip susceptibility.This requires the use of a coupled thermo-hydromechanical (THM) model to simulate the subsurface fluid injection process and associated stress changes, in particular the magnitude of perturbing normal and shear stresses and their orientations relative to the vulnerable fracture planes of investigation.The authors (Cao et al. 2022) developed such a fully coupled THM model for the Húsmúli area to examine the relative contributions of direct pore pressure, poroelastic stressing, and thermal stressing to the potential for seismicity, and such a model could be used for probabilistic evaluation of seismicity due to various causal mechanisms.
The seismic occurrence probability forecasted is influenced by both the spatial distribution of fluid pressure change and the density of underlying fractures (Eq.9).In this work, the fluid circulation at the Hellisheiði geothermal site has reached steady-state flow conditions during the periods of investigation, therefore, the estimated probability for induced seismicity is predominantly influenced by the varying abundance of potentially vulnerable underlying fractures under the influence of prevailing stress and fluid pressure over the fluid injection period.Nevertheless, this probabilistic seismic evaluation method is also applicable to scenarios with rapid-changing fluid pressures (e.g., increasing injection rate, cyclic injection operations, and transient shut-in), provided that the spatiotemporal evolution of fluid pressure change can be faithfully captured by the hydrogeological modelling.In this regard, the methodology proposed can be used to evaluate future fluid injection scenarios and to identify optimal operation conditions in terms of the location, duration and rate of fluid injection.
It is noteworthy that the seismic occurrence probability estimated is greatly dependent on the considered scale in time and space (Langenbruch et al. 2018).A larger gridblock dimension increases the forecasted seismic occurrence probability.Likewise, a longer timescale leads to a larger forecasted seismic occurrence probability within a certain region.The forecasting period selected should neither be too large to ensure the timeliness of seismic forecasting, nor too small to avoid low data statistics owing to the limited number of seismic catalogues (Cao et al. 2020a).In this work, the forecasting period selected is half a year, comparable to that used in seismic evaluation for other fluid injection field sites (e.g., 1 year in Langenbruch et al. 2018).

Conclusions
Fracture criticality, defined as the gradient of critical fluid pressure change to trigger seismicity (Δp c /h), was proposed as a novel metric of fault slip susceptibility.The fracture criticality is an appropriate metric to characterise perturbations required to trigger seismicity (or the critical state of pre-existing fractures) in 3D space, as its variability is independent of reservoir depth, and only reflects uncertainties involved in fracture attributes and rock properties.The fracture criticality at the Hellisheiði geothermal site was interpreted using integrated analysis of seismic observations and hydrogeological simulation of fluid injection operations in the faulted reservoir at this site.It has been found that the fracture criticality within both fault and off-fault zones shows natural variability (mostly ranging between 0.001 and 2.0 bar/km), and the values estimated roughly follow Gaussian distributions.Fault zones tend to be characterised by larger fracture criticality values than off-fault zones.Regardless of the spatial location of the region, the range of fracture criticality values within the region varies distinctively for different fluid re-injection periods.
Based on statistics of the fracture criticality, a probabilistic seismic susceptibility evaluation framework was developed to evaluate the potential for fluid injection-induced seismicity in both faulted and off-fault reservoirs by combining seismic observations and hydrogeological simulation.The framework considers the fluid pressure increase driven by fluid injection operations, the intrinsic variability of fracture criticality, and regional fracture density related to neighbouring fault structures.This methodology was applied to evaluate susceptibility to induced seismicity at the Hellisheiði site, and relatively high probability of seismic event occurrence was estimated for fault zones around the five re-injection wells.Areas of high estimated fracture slip susceptibility experienced high levels of induced seismicity over the fluid re-injection period of investigation.During the 2-year re-injection period investigated, whilst F1, F2 and F3 fault zones evolved seismotectonically from an active state during the first two half-year periods to a much suppressed state during the last two half-year periods, the seismotectonic state of F4 and F5 fault zones and off-fault zones maintained fairly low levels throughout.The seismotectonic state estimated for each region during each half-year modelling period is generally consistent with the seismic susceptibility forecasted based on the statistics of fracture criticality.

Fig. 1
Fig. 1 Schematic illustration of critical fluid pressure changes Δp c and Coulomb failure stress changes ΔCFS c required to cause the slippage of a fracture plane at different reservoir depths.The outward unit normal vector of the fracture plane makes an angle of β with the direction of the maximum principal stress

Fig. 2
Fig.2Procedure for probabilistic fluid injection-induced seismicity forecasting using integrated monitored seismicity interpretation and hydrogeological modelling

Fig. 3 a
Fig. 3 a Map of the Hengill volcano and the Hellisheiði and Nesjavellir geothermal fields showing elevation contours, surface faults/ fractures, and main production and re-injection wells and well paths projected to the surface (after Tómasdóttir 2018).b Map of the The reservoir model constructed was based on a geological model of the Hellisheiði geothermal field, which incorporated surface mapping data, well data, well logging data, geophysics data, and laboratory analyses data (Gunnarsdóttir and Poux 2016).The geological model is comprised of a combination of a lithological model and a structural model.The reservoir model developed consists of two main types of geological formations, the hyaloclastic formation (mostly at 0 ~ 1600 m depth) and basalt basement (mostly at ~ 1600-3000 m depth).

Figure 8
Figure 8 presents the histograms of critical fluid pressure change and gradient of critical fluid pressure change

Fig. 6
Fig. 6 Spatial distribution of fluid pressure in the reservoir model of the Húsmúli area.Coordinates of the model geometry follow the ISN93 geographic coordinate system suitable for use in Iceland

Fig. 7
Fig. 7 Geothermal fluid injection rates and history matching of field bottomhole pressure time series at 5 re-injection wells at the Húsmúli area during 15 April-15 October 2020

Fig. 8
Fig. 8 Histograms of critical fluid pressure change Δp c and gradient of critical fluid pressure change Δp c /h (fracture criticality) for induced seismicity recorded at different areas at Húsmúli during 15

Fig. 9
Fig. 9 Susceptibility to fluid injection-induced seismicity estimated based on statistics of fracture criticality at 1 km depth at the Húsmúli area for the next half-year fluid re-injection period (16 October 2020-15 April 2021): a spatial distribution of modelled re-injection

Fig. 12
Fig. 12 Evolution of the seismotectonic state of different regions (F1-F5 and off-fault zones) at Húsmúli: a the density of seismicity against the average fracture criticality, b the seismotectonic state over time, and c the seismogenic index over time.Four points connected

Table 1
Mohr-Coulomb fracture slippage criterion-based fracture slip susceptibility indices *τ is the shear stress along the fracture, and σ and p are the normal stress and fluid pressure applied on the fracture plane.f is the friction coefficient of the fracture, and α is Biot's coefficient (α = 1 in this work) or 3D fracture network representation (Tómasdóttire 2018) to match field tracer test results.The latter is considered more suitable in that 3D fracture network representation could represent fluid flow in different directions and better account for hydrological effects imposed by far-field injection and production wells.Recently, based on outcrop mapped fracture systems and tracer test results at the Hellisheiði geothermal field,Mahzari et al. (