Pseudo‐Prospective Forecasting of Induced and Natural Seismicity in the Hengill Geothermal Field

The Hengill geothermal field, located in southwest Iceland, is host to the Hellisheiði power plant, with its 40+ production wells and 17 reinjection wells. Located in a tectonically active area, the field experiences both natural and induced seismicity linked to the power plant operations. To better manage the risk posed by this seismicity, the development of robust and informative forecasting models is paramount. In this study, we compare the forecasting performance of a model developed for fluid‐induced seismicity (the Seismogenic Index model) and a class of well‐established statistical models (Epidemic‐Type Aftershock Sequence). The pseudo‐prospective experiment is set up with 14 months of initial calibration and daily forecasts for a year. In the timeframe of this experiment, a dense broadband network was in place in Hengill, allowing us to rely on a high quality relocated seismic catalog. The seismicity in the geothermal field is characterized by four main clusters, associated with the two reinjection areas, one production area, and an area with surface geothermal manifestations but where no operations are taking place. We show that the models are generally well suited to forecast induced seismicity, despite some limitations, and that a hybrid ETAS model accounting for fluid forcing has some potential in complex regions with natural and fluid‐induced seismicity.


Introduction
Iceland's capital region gets most of its district heating and a third of its electricity from the exploitation of neighboring geothermal plants (Gunnlaugsson & Ívarsson, 2010).Located 30 km east of Reykjavik, the Hengill field hosts two major power plants supplying the capital region: Hellisheiði and Nesjavellir, with installed capacities of 303 MW e plus 133 MW th and 120 MW e plus 150 MW th respectively (Hersir et al., 2009).Nesjavellir was commissioned in 1990 while Hellisheiði started production in 2006, with the subfield of Hverahlíð beginning full scale production in late 2017.
Iceland is a seismically active country, especially along the mid-Atlantic ridge that crosses through the island (Jakobsdóttir, 2008).On top of this natural seismicity, the exploitation of deep geothermal energy can cause induced seismicity (Grigoli et al., 2017).In the Hengill field, numerous episodes of seismicity have been recorded since the early 1990s with the start of instrumental catalogs in southwest Iceland (SIL network, Icelandic Met Office).Both volcano-tectonic and induced sequences illuminating the field have been mapped to shallow depth, and the discrimination between natural and induced seismicity is particularly difficult in this area.
Although relatively isolated and with a low building density, the area around the Hengill field has seen several large events in the last decades, including two likely induced events of magnitude around 4 in October 2011.The operators of the geothermal plants thus need tools to assess and mitigate the seismic risk posed by their injection and production operations.In this study, we use probabilistic models to forecast seismicity happening in the Hellisheiði field between late 2018 and early 2021.
Statistical models have shown some encouraging results to model induced seismicity (Király-Proag et al., 2016;Mancini et al., 2021;Verdon & Budge, 2018), although these models do not account for coupled processes in the subsurface.We use two different classes of statistical models developed for very different use cases: Seismogenic Index, developed for injection-induced seismicity; and Epidemic-Type Aftershock Sequence (ETAS) models that are widely used for natural seismicity modeling.The Seismogenic Index is a simple yet robust model that relies on a linear relationship between volume rate and seismicity rate, and has proven its reliability in numerous induced seismicity sequences (Mignan et al., 2017(Mignan et al., , 2021;;Shapiro et al., 2010).Recent work showed that its parameters can be updated in near real time during stimulations (Broccardo et al., 2017(Broccardo et al., , 2019;;Mignan et al., 2019).The Epidemic-Type Aftershock Sequence model relies on several empirical observations, from the Gutenberg-Richter power law, to Omori's law for aftershock decay.ETAS models are at the forefront of operational earthquake forecasting (Marzocchi et al., 2014;Nandan, Kamer, et al., 2021).
We try to answer a few questions: How do Seismogenic Index-type models (tailored to induced seismicity) compare to ETAS based models (adapted from natural to mixed seismicity with fluid forcing term)?How do these models perform in forecasting the seismicity, retrospectively?What can these models tell us about the seismicity in the Hengill geothermal field?

Geological Context
The Hengill geothermal field is located in southwest Iceland, 30 km east of Reykjavik.The field is nested between three volcanoes: Hengill to the north, Hrómundartindur to the north-east, and Graensdalur to the east.Mount Hengill was last active around 2000 years ago, while Hrómundartindur was last active circa 10,000 years ago, and Graensdalur has been extinct for 300,000 years (Foulger & Toomey, 1989;Jousset et al., 2011;Sánchez-Pastor et al., 2021).The area forms the junction of three tectonic systems: The Reykjanes Peninsula oblique rifting system (RP), the South Iceland Seismic Zone (SISZ) and the West Volcanic Zone (WVZ) (Tomasdóttir, 2018; Figure 1 lower left panel).
Located to the south of Mount Hengill, the Hellisheiði subfield is characterized by a water dominated fractured system, comprising basaltic lava layers, hyaloclastites series and dyke intrusions (Franzson et al., 2005;Snaebjörnsdóttir et al., 2018), with and average porosity of 10% (Gunnarsson et al., 2011).Only the southernmost area of Hverahlíð presents a different composition with mostly lava series.The formation temperature averages between 220 and 250°C at 1000 m b.s.l.(Gunnarsson, 2013).The produced fluid is a water-steam mix ranging in temperature between 240 and 320°C.
The area is characterized by the Hengill fissure swarm, which strikes 30°N forming a graben structure around 40 km in length (Saemundsson, 1992).Extensional structures (normal faults) and eruptive fissures are a common occurrence (Steigerwald et al., 2020), as well as strike-slip faults oriented N-S associated with the SISZ.
The field has seen a major volcanic uplift event between 1994 and 1999, with up to 8 cm of inflation, close to 100,000 earthquakes recorded with a largest magnitude of M L 5.5 (Blanck et al., 2021;Jakobsdóttir, 2008).Two earthquakes of magnitude 6 have been recorded in 2008 in the eastern part of the field, reactivating strikeslip faults (Decriem et al., 2010).Following the late 1990s uplift event, subsidence has been observed since the mid 2000s (Ducrocq et al., 2021), coinciding roughly with the commission of the Hellisheiði power plant in 2006.

Areas of Interest
The Hengill geothermal field is host to two large power plants: Nesjavellir in the north, and Hellisheiði in the south (Figure 1).In this study, we only focus on the southern part of the field, which comprises of various production areas as well as two deep reinjection areas.We selected four areas of seismological interest, based on the spatial clustering in the COSEISMIQ catalog.These areas are highlighted by the colored boxes in Figure 2.

Production Areas
The Hengill field produces a water-steam mixture at an average rate of 38Mton/year (for the 2012-2015 period, Juncu et al. (2017)).The production occurs in 40+ wells located in two distinct areas: The central Hellisheiði field, and Hverahlíð, a sub-field to the south.

Hellisheiði:
The main production area of the Hellisheiði field has not been associated with much seismicity since the end of the volcanic event in 1999.The drilling of the production wells since 2001 in the area did not trigger seismicity (except for well HE-08).However, since 2009 and the end of drilling, there has been sustained levels of seismicity, with some seismic swarms happening intermittently (Hjörleifsdóttir et al., 2021).The area is however rather quiet during our period of interest and is thus not looked at in detail in this study.
Hverahlíð: Hverahlíð is located in the south of the Hengill region and is considered a subsystem of Hellisheiði.It hosts the most powerful wells of the Hengill field and is associated with a shallow (2-3 km depth) low velocity anomaly (Sánchez-Pastor et al., 2021).The area is also characterized by a relative abundance of strike-slip faults linked to the SISZ, making the area likely even more permeable than the rest of the field (Franzson et al., 2010).The drilling of the production wells in Hverahlíð was not associated with seismicity except for well HE-21 in 2006.The production in the area started in fall 2017, coinciding with an increase in the seismicity which remains elevated as the production continues (Hjörleifsdóttir et al., 2021).

Reinjection Areas
The reinjection of spent fluids in a geothermal field is mandated for multiple purposes: To sustain reservoir pressures, avoid subsidence and reservoir compaction, enhance the natural recharge of the system, improve thermal extractions along flow-paths, and comply with environmental regulations (Axelsson, 2012).In Hellisheiði, the reinjection occur mostly in two clusters of wells (Húsmúli and Gráuhnúkar) with a few other reinjection wells scattered within the production areas.Gráuhnúkar: Gráuhnúkar is the first dedicated reinjection area of the Hellisheiði field, and is located in the southwest of the region.The six injection wells were commissioned between 2006 and 2008, and are continuously used for reinjection purposes since then (Hardarson et al., 2010).The formation temperature is much higher than expected with temperatures reaching up to around 300°C.The start of reinjection in Gráuhnúkar did not lead to an alarming increase in the seismicity rate, which remained low until the injection rate was increased in mid-2011.Since 2011, the seismicity has been low to moderate and seems to correlate with the injection rates and to take off when the 300 L/s mark is reached (Flóvenz et al., 2015;Ritz et al., 2021).
Húsmúli: The drilling of the eight wells in Húsmúli started in 2007 and continued until mid-2011.The reinjection started systematically in five of the wells in September 2011 (although some injection had been done in one of the wells since 2009 (Gunnarsson et al., 2015;Hardarson et al., 2010)).Seismicity started during the drilling operations, linked to repeated circulation losses, and remained very high during the first nine months of reinjection until the injection rates were reduced (Ágústsson et al., 2015;Gunnarsson, 2013;Hardarson et al., 2010;Kristjansdottir et al., 2021;Ritz et al., 2021).During the first year of reinjection, a 2 cm uplift was observed in Húsmúli (Juncu et al., 2018).A shallow (∼3 km) deflating source has been observed in the late 2010s, potentially linked to the localized fluid extraction and circulation (Ducrocq et al., 2021).

Ölkelduháls Area
Ölkelduháls is located to the east of our region of interest.It is characterized by surface geothermal manifestations like hot springs and fumaroles probably linked to the residual heat of Hrómundartindur volcano.A high velocity anomaly at depths 2-4 km has been identified (Jousset et al., 2010).In 1995, an exploratory well was drilled in the area, and two more have been drilled since but are not active.Ölkelduháls is very close to the center of the late 1990s volcanic uplift (estimated between 5 and 6 km depth), experienced subsidence between 2006 and 2017 and was at the center of the an inflating source mapped in 2017-2018 with up to 12 mm vertical displacement (Ducrocq et al., 2021).These deformation cycles are not well understood but have been theorized to be the results of either irregular magmatic intrusions or hydrothermal fluid migrations (or a combination of both with different layers seeing different fluids intrude or migrate).Ölkelduháls is the only region that is not highlighted because of its suspected induced seismicity, but quite some seismicity is recorded in the area.

Seismic Catalog
The seismic data was acquired as part of the COSEISMIQ project (http://www.coseismiq.ethz.ch/en/home/)during which the Hengill region was instrumented with a small aperture seismic array between December 2018 and February 2021.Multiple catalogs were created with different quality thresholds and relocation techniques.
Here we use the high quality catalog with relative relocation (Grigoli et al., 2022).This catalog contains about 8,500 events distributed in a 35 × 30 km area centered around the Hellisheiði power plant (Figure 1).We use all the events in the catalog to determine the statistical properties shown in Figure 3.The active area highlighted in the main panel of Figure 1 corresponds to the limits of Figure 2.
The b-value and magnitude of completeness of the catalog are estimated jointly using a method proposed in the literature by Clauset et al. (2009) and Mizrahi et al. (2021a).We assume the magnitude of completeness (M c ) to be constant both in time and space in the area of interest and test different values of M c , calculating for each M c the corresponding b using the maximum likelihood approach (Marzocchi & Sandri, 2003).We then check whether the observed cumulative magnitude distribution function is plausible to be a realization of the fitted Gutenberg-Richter (GR) law as follows.The cumulative distribution function (CDF) corresponding to the fitted discretized power law is compared to the observed cumulative distribution function using the Kolmogorov-Smirnov distance (KS-distance).The probability p M c of observing a KS-distance of at least D M c from a sample randomly drawn from a discretized GR power law is estimated through the simulation of 10,000 of such random draws.p M c is defined as the fraction of simulated samples for which the KS-distance is at least D M c .The smallest M c for which this probability is greater than or equal to 0.1 is considered to be the magnitude of completeness of the catalog.From this analysis we conclude that the magnitude of completeness of the catalog is 0.3 (Figure 3a), and that the b-value is equal to 0.93 ± 0.01 (Figure 3b).

Hydraulic Data
The Hellisheiði geothermal field counts 17 active reinjection wells and around 45 active production wells at all times (Figure 2).These wells are directional wells for the most part, but we only have access to the wellhead coordinates for the production wells, which we used as proxy for the feedzone locations.
The Hengill field has been in continuous operation since 2006, and with more than a decade of continuous circulation, one can assume that the volume of the reservoir is at its maximum extent by the start of the COSEISMIQ project.For this reason, the injection and production rates per well obtained from the power plant operator (Reykjavik Energy) are distributed in a symmetrical bi-variate Gaussian distribution around the wellhead to account for fluid migration in the subsurface with daily granularity.The covariance matrix of the bivariate Gaussian is defined as Σ = 5 × 10 5 0 0 5 × 10 5 ( ) (in degrees longitude/latitude).This value of 10 5 chosen for the variance stems from a comparison of the extent of the reservoir using a regional TOUGH2 model (Bjornsson et al., 2003;Gunnarsson et al., 2021), and has been slightly increased to account for the use of the wellhead as a proxy for the feedzones of the know-to-be deviated wells.The produced and injected volumes are summed up into a compound volume that we will use as input for the models taking in fluid volume or rate.
It is worth noting that this isotropic distribution of the volumes is only a coarse proxy and doesn't exactly match with the areas of seismological interest (Figure 4).In particular, the volumes distributed in Gráuhnúkar are offset to the north-west compared to the seismic cluster.The wells in Gráuhnúkar actually strike southeast, and line up pretty well with the seismicity cluster highlighted by the box in Figure 4 (Hjörleifsdóttir et al., 2021;Ritz et al., 2021).More physical approaches to the distribution of fluids are rendered complicated by the two-phase nature of fluids in the reservoir.Indeed, one cannot simply rely on the hydraulic diffusivity as the high portion of gas-phase would result in a highly unrealistic distribution of the fluids.

Models
Two families of models are implemented: a Seismogenic Index type and a class of Epidemic-Type Aftershock Sequence models.

SI-Seismogenic Index Model
This model is based on a purely empirical relationship between the injected/produced volume and the seismicity, and contains a set of parameters describing the geological and seismological characteristics of the targeted site.
Various versions exist in the literature after the pioneering work of S. Shapiro's group (Shapiro, 2018).Other forms account, for example, for the injection rate rather than the injected volume, as well as for an exponential decay of the seismicity after shut-in (Broccardo et al., 2017;Mignan et al., 2017).The proposed model follows the approach by Broccardo et al. (2017).The expected number of events N expected follow the relationship: where V inj/prod is the compound (i.e., injected and produced) volume, a f b the activation feedback (in m 3 ), and b the b-value of the Gutenberg-Richter power law.This underground activation feedback a f b is the a-value of the Gutenberg-Richter power law normalized by the volume such that a f b = a log 10 (V), and is also known as the seismogenic index in the poroelastic context (often noted Σ; Shapiro et al. (2010); Dinske and Shapiro (2013)).This parameter has also been interpreted as resulting from geometric operations on a static stress field produced by a change in the volume in the underground (Mignan, 2016).For simplicity, we refer to the model as the 'Seismogenic Index' and the productivity parameter a f b as the activation feedback, and do not assign any specific physical meaning beyond that of fluid-earthquake productivity factor.
Equation 1 can be translated into an equation to determine the seismicity rate for active phases as: In this model, the linear relationship between V (t) and λ(t,m > M c ) derives from the linear relationship between injection rate and overpressure, and neglects potential temporal changes in injectivity (Mignan, 2016).Note that this equation has been extensively tested on injection-induced sequences, but not on the production phase (Broccardo et al., 2019;Mignan et al., 2017), although the classical Seismogenic Index approach has been extended to production cases (Shapiro, 2018).
We do not account for post-injection phases and the associated seismicity decay-as modeled by Mignan et al. (2017)-as the Hengill field is in constant operation, both for production and injection.

Parameter Calibration
For the estimation of the parameters of the Seismogenic Index model, M c is taken from the catalog, while a f b and b are optimized for each cell by minimizing the log-likelihood function (Equation 3).The Maximum Likelihood Estimate is modified from Broccardo et al. (2017) to only account for active injection and/or production phases.
ln This data-driven cell by cell approach allows the model to assume that the underground feedback parameter a f b varies spatially within the area of interest.For cells with less than 50 events above completeness, the b-value of the catalog (b = 0.93) gets assigned, and the underground feedback parameter a f b gets calibrated using the minimization of the log-likelihood function.The variability in a f b and b is small from one cell to the next, but significant between different subregions, which indicates different stressing and/or fluid-phase conditions in the field (see Figure S1.1 in Supporting Information S1 for maps of the optimized a f b and b parameters at the end of the learning period; note the uniform b-value in the central region denoting that less than 50 events are counted in the cell at the end of the learning period).
Figure S1.2 in Supporting Information S1 gives an impression of the ability of the model to fit to the data by showing the difference between the recorded and expected events during the learning period.The Seismogenic Index model fits very well in the Húsmúli, Hellisheiði, and Hverahlíð regions, where there are both volumes injected and/or produced and seismic events.The model is not able to fit anything in Ölkelduháls as there are no injection or production activities.The fit in Gráuhnúkar is good in the eastern part where the volumes are present, but unsatisfactory to the west where there are no volumes distributed.
In this current state, the Seismogenic Index model is probabilistic by nature and carries the assumption that induced seismicity can be described as a Poissonian process.From the optimized parameters, we calculate the mean rate of seismicity in the cell at this given time and extrapolate a Poissonian distribution from it to be able to make a probabilistic comparison of the forecast to the recorded seismicity.

ETAS-Epidemic-Type Aftershock Sequence Models
ETAS models (Ogata, 1988) view seismicity as a combination of background and triggered earthquakes.Background earthquakes result from tectonic forces or anthropogenic factors such as fluid injection and/or production.These background earthquakes can trigger a cohort of aftershocks which can then trigger their own aftershocks and so forth.In its simplest form the ETAS model describes the conditional seismicity rate of magnitude m events, λ t,x,y,m|H t ( ), at any location (x,y) and time t as where β = ln(10) • b is the GR-law exponent in base e, μ is the background intensity function, which may or may not depend on space or time, H t is the history of the process up to time t, and g(m,Δt,Δx,Δy) describes the rate of aftershocks triggered by an event of magnitude m, at a time delay of Δt and a spatial distance (Δx, Δy) from the triggering event.Here we use the specific kernel g defined as: see Nandan, Kamer, et al. (2021) and Mizrahi et al. (2021b).
This kernel combines three components to describe the aftershock rate of an event of magnitude m at a temporal and spatial distance of (Δt, Δx, Δy) from its origin.The productivity law which quantifies the number of events directly triggered by an event of magnitude m i as k 0 e a m i M c ( ) , where a parameterizes the magnitude dependent part of the aftershock productivity, and k 0 parameterizes the magnitude independent part.The time-based Omori-Utsu aftershock decay is given as e Δt/τ (Δt+c) 1+ω , where c avoids a singularity at time 0 and has been brought in relation with different phenomena ranging from short-term aftershock incompleteness to physical mechanisms involved in earthquake nucleation, 1 + ω quantifies how fast the aftershock rate decays with time, and τ determines the onset of the exponential taper that allows negative ω values to be obtained.The spatial aftershock decay is described through an isotropic power law, Δx 2 + Δy 2 ) 1+ρ) , where 1 + ρ describes how fast aftershock rate decays with distance from the triggering event, and d ⋅ e γ m M c ( ) describes a magnitude dependent distance at which the aftershock rate decay starts.Note that it is implicitly assumed that the completeness magnitude M c coincides with the magnitude of the smallest event which is capable of triggering other events.

Calibration of Parameters
To calibrate the parameters of the ETAS model, we apply the Expectation-Maximization algorithm (Veen & Schoenberg, 2008).In this iterative algorithm, one starts with a random initial guess for each of the parameters that need to be calibrated.The expectation step and the maximization step are repeated, updating the current estimation of the parameters, until the difference between the parameters of two iterations falls below a desired threshold.In the expectation step, the probabilities p ij that event e j was triggered by event e i , the probability p ind j that event e j is independent, the expected number of background events n, and the expected number of directly triggered aftershocks li of each event e i , given the ETAS parameters of the current iteration, are estimated as and Here, g kj = g(m k , t j t k , x j x k , y j y k ) is the rate of aftershocks of event e k at the location and time of event e j .With these definitions, the independence and triggering probabilities are proportional to the contribution of background and aftershock triggering terms at any given time and location.
In the Maximization step, the ETAS parameters are optimized to maximize the complete data log-likelihood, and these optimized parameters are used in the next iteration of the Expectation step.

Issuing a Forecast Through Simulations
Once the ETAS parameters are calibrated based on the training catalog, a forecast is issued by simulating 100,000 possible continuations of this training catalog.This includes the simulation of aftershock cascades of the events in the training catalog, and the simulation of background earthquakes which fall into the forecasting period plus their cascades of aftershocks.A detailed description of the simulation algorithm is given in Mizrahi et al. (2021a).
Training as well as simulation are done for the full region of the catalog, although they are evaluated only in the active geothermal area (Figure 1).This ensures that aftershock triggering which goes beyond the borders of the relatively small active area is still captured by the model.
Based on the simulated catalogs, the likelihood p(k) of k events to occur in a spatial cell of interest during the forecasting period is then given through the empirical distribution as where n(k) is the number of simulations for which k events are observed, m 0 is the number of values k between zero and 100 for which n(k) = 0. To avoid zero probabilities for values of k that do not appear in the simulations, that is, n(k) = 0, this definition of p(k) includes a water-level probability for all k ≤ 100.

ETAS Variants
The four ETAS variants applied in this study differ only in their formulation of the background seismicity term.

ETAS-0
In the most basic variant (standard ETAS), ETAS-0, μ is considered constant in space and time during parameter calibration.This means that the background term which determines independence and triggering probabilities in Equations 6 and 7 is the same for all events, irrelevant of where and when they take place.
For the simulation of catalog continuations, the background seismicity term is space-dependent.While the number of simulated background earthquakes during the forecasting period is simulated from the constant parameter μ, their locations are simulated by randomly drawing locations of earthquakes in the training catalog, weighted by their probability of being a background event, and adding a distortion drawn from a Gaussian distribution with mean 0 and σ = 0.5 km.

Varying Background Rate (ETAS-v)
The second ETAS variant considered in this study uses a space-varying background term μ(x,y) during model calibration and is designated afterward as ETAS-v.In each iteration of the Expectation-Maximization algorithm, μ(x,y) is calculated as where k(Δx j , Δy j ) is a Gaussian kernel with bandwidth σ = 0.5 km applied to the distance (Δx j , Δy j ) of event e j to the location (x,y), and T is the time length of the training catalog.This model corresponds to the flexible ETAS model with free background described by Mizrahi et al. (2023).

ETAS With Fluid Forcing (ETAS-f)
ETAS-type models have been used in induced seismicity contexts (Bourne & Oates, 2017;Mena et al., 2013), generally in regions with low natural background seismicity where the calibration could focus on induced earthquakes only.However, fluid-driven seismicity has distinct spatio-temporal characteristics from tectonicloading driven seismicity and requires its own parameters when we want to model complex areas with both induced and natural seismicity.
For single-well injections, Bachmann et al. (2011) introduced an external forcing on the background rate linearly proportional to the injection rate, and found this model to perform better than the standard ETAS in a pseudoforecasting experiment of the hydraulic stimulation in Basel (Switzerland).More recently, the same approach has been applied to a hydraulic fracturing context with promising results (Mancini et al., 2021).
We here introduce an ETAS variant ETAS-f with the background term which, in addition to the tectonic background rate μ tect , comprises a fluid background rate μ f luid = ι • V(x,y,t) which is proportional to the compound volume (sum of produced and injected volumes) V(x,y,t).The independence (or background) probability p ind j defined in Equation 7 can be split into a tectonic and a fluid part and analogously the expected number of background events is the sum of the expected number of tectonically triggered and fluid triggered background events, ntect and nfluid .The new parameter ι events m 3 [ ] can be calibrated in the Maximization step of the Expectation-Maximization algorithm as When a forecast is issued for this ETAS-f model, two types of background earthquakes are simulated, tectonic and fluid-induced ones.Both types of background events can trigger cascades of aftershocks (the aftershock triggering parameters are the same for both tectonic and fluid-induced events).The number of induced earthquakes simulated in a given cell on a given day is determined by the compound volume of that cell on that day.Note that this is based on the assumption that the planned compound volume for that day is already known at the time the forecast is issued, which is a reasonable assumption for short enough forecasting horizons and in a geothermal field with stable exploitation conditions like Hellisheiði.The locations of simulated induced earthquakes within a cell are generated using a uniform spatial distribution inside the cell for which the volume is given.
This modification of the ETAS model can be applied to the standard ETAS version as well as the version with spatially varying (tectonic) background seismicity.This yields a total of four variants of ETAS: ETAS-0, ETASv, ETAS-f, and ETAS-fv.
Note that this formulation of ETAS with fluid forcing is based on several simplifying assumptions.No difference is made between injected and produced volume, the effect of fluids on seismicity is assumed to be immediate and only valid on the current forecasting horizon, induced seismicity is assumed to be proportional to compound volume, etc.If this simple approach produces promising results, these assumptions can be revisited in future studies and a more realistic formulation of the fluid-induced term can be applied.In particular, production and injection could be split into two terms to better account for the differences in physics behind the triggering of induced events in the different settings.

Model Comparison Framework
Evaluating the performance of a model to reproduce or forecast induced seismicity is not an easy task, but is a necessary step to the implementation of ensemble modeling with model-specific weights (Király-Proag et al., 2018).Different classes of models have intrinsic assumptions that make it difficult to do direct comparison, in particular in the statistical or stochastic expression embedded in the models.Guidelines have been proposed to standardize the evaluation of models, in particular the "Induced Seismicity Testbench" (Király-Proag et al., 2016), and the guidelines developed by the Collaboratory for the Study of Earthquake Predictability (CSEP, https://cseptesting.org/).
Each model provides a forecast consisting of a probability distribution of the number of earthquakes for each time and spatial bin.To evaluate and rank model performance, we use a probabilistic score: The log-likelihood, which is the natural logarithm of the probability of the model matching the recorded number of events in the time and space bin.The Seismogenic Index model doesn't yield a forecast in "non-fluid" cells, so it is only evaluated in a portion of the cells (colored cells in Figure S1.1 in Supporting Information S1).This practically means that a score of "0" is assigned to cells without a forecast, thus allowing us to compare on the entire region even though not all models yield a forecast for each cell.This approach could be criticized as it assumes what is effectively a perfect forecast in cells where no forecast is produced, but the Seismogenic Index model is only designed to model fluidinduced seismicity and it would be unfair to penalize it in non-fluid cells.For each time bin, the log-likelihood is summed up over the spatial cells to obtain the "score" of the model at this forecasting horizon.The log-likelihood is then accumulated over time into the cumulative spatial joint log-likelihood to show the evolution of the forecasting performance and compare models.
We use the Information Gain to measure the predictive performance of a model compared to that of a reference model (Kagan & Knopoff, 1987).In this study, the standard version of ETAS is used as the reference model.The Information Gain is calculated as the difference of the log-likelihoods of the models, and is positive when the model performs better than the reference model and vice-versa.
ETAS models provide a full probability distribution of the expected event number per spatial cell as described in Equation 10.In the case of the Seismogenic Index, we assume a Poissonian distribution with the forecast number of event as the mean λ of the distribution.

Forecasting Framework
The pseudo-prospective forecasting experiment runs between December 2018 and the end of January 2021, the period being covered by the COSEISMIQ high-quality relocated catalog.The first 14 months are used as a learning period to train the models (1.12.20218-31.1.2020).We then perform daily forecasts for 366 days (1.2.2020-31.1.2021).After each forecasting horizon, the day is added to the calibration data set (Figure 5).
The area of interest is centered around the production region of Hellisheiði and encompasses the reinjections regions of Húsmúli and Gráuhnúkar, as well as the production sub-field of Hverahlíð.We discretize the space in 0.005°latitude × 0.005°longitude spatial bins (rectangular cells measuring approximately 560 × 240 m; shown in the background of Figure 2).The forecasts are issued for each spatial and temporal bins.Subsets of the forecast are examined for the regions of Húsmúli, Gráuhnúkar, Hverahlíð, and Ölkelduháls (boundaries defined in Figure 2).

General Performance of the Models
We evaluate the absolute and relative predictive performance of the models over a one year period.Figure 6 shows a comparison of the cumulative number of events recorded (dashed line) and forecasted by the different models (color-coded), for the active region (top panel) and for the subregions of special interest.The Seismogenic Index model has a very linear trend, resulting from its low number of parameters, whereas ETAS models are more flexible and display variations in the modeled seismicity rate.Overall, ETAS-0 and ETAS-v show a closer fit to the total number of events, even though the seismicity rate is overestimated between March and late October 2020.The ETAS models with fluid-forcing and the Seismogenic Index have a better fit initially, but fall behind in the second half of the forecasting experiment period.This discrepancy between ETAS-0 and ETAS-f models is a consequence of the lower branching ratio η in the ETAS-f versions, which means that less events are interpreted as aftershock.This is a direct consequence of the doubling of the background terms (tectonic μ tect and fluid μ f luid ) that leads to more flexibility of the model to accommodate events as background (tectonic or fluid).Thus, even though the background rate is high, fewer events are simulated because aftershocks are missing compared to the non-fluid versions.
Figure 7 shows the cumulative spatial log-likelihood in time for each model for the active region.This metric evaluates the consistency of the models with the observations, with a log-likelihood closer to zero indicating a better model performance.All models exhibit a somehow jagged behavior during periods of higher seismic intensity, which is expected.The Seismogenic Index is quite competitive with the ETAS models until October 2020, when it plunges down and is unable to explain a sudden increase in seismicity, that the ETAS models manage to accommodate.The cumulative number of events presented in Figure 7 only shows the events recorded during the forecasting period, to account for all the recorded seismicity, and does not include the 1,093 events above completeness recorded during the learning period.
In Figure 8, we compare the different ETAS variants and the Seismogenic Index to a reference model, by visualizing their cumulative information gain over time with respect to the standard version ETAS-0.It is worth noting that the vertical scale of the figure is cut-off at 600, but that the Seismogenic Index's information gain plunges all the way down to 2900 by the end of the forecasting period.The Seismogenic Index performs quite well in the beginning, outperforming all ETAS models.Its downfall comes in October 2020, as the first highintensity period takes place with sudden high rates of seismicity that the model cannot accommodate.ETAS-f outperforms the standard ETAS model by the end of the forecasting experiment, showing that in a geothermal area with active injection and production, accounting for the volumes does improve the forecasts.However, the gain is limited for most of the study period until October 2020 which sees a shift in the baseline seismicity rate with heightened activity in the reinjection area of Húsmúli.On the other hand, the varying background versions of ETAS (ETAS-v and ETAS-fv) don't perform as well as the standard version of ETAS.Possibly, the formulation of the varying background seismicity used here is not well-suited to describe a complex field like the present one.
Giving too much flexibility to the background part of seismicity can make the model wrongly interpret local phenomena or fluid-induced seismicity as variations in background seismicity rate.In future studies, the bandwidth σ determining the spatial smoothing in Equation 11 could be calibrated as a hyper-parameter, or a different type of smoothing could be tested to obtain better performance.
When comparing Figures 6-8, one might be puzzled as it appears that models getting the closest fit in terms of total number of events are not ranked higher in the log-likelihood evaluation.However, the model matching the forecasted cumulative number of events the closest is not guaranteed to have the best performance.Indeed, by estimating the log-likelihood and accumulating it over time, we jointly evaluate the ability of a model to fit changes in the seismicity rate.In Hverahlíð, the Seismogenic Index is closest to reality in terms of cumulative forecasted number of events (Figure 6 bottom right panel), but is performing worse than all ETAS models (Figure 7 bottom right panel), which is reflected as an increasingly more negative information gain (Figure 8 bottom right panel).This is the result of ETAS models being better at reproducing temporal variations in the seismicity rate.

How Do Models Perform in the Subregions of Interest?
We evaluate the models in four separate subregions of specific seismological interest by summing the loglikelihoods within the boundaries of each region.This evaluation for subregions is done from the outputs of the entire field for each model that is, the models are not trained specifically on the subregions.
Ölkelduháls is not a region with active injection or production so the Seismogenic Index Model does not give any forecast for the region (Figure 7 middle right panel).Seventy-three events were recorded during the learning period, and 69 events are recorded during the forecasting period.In this subregion, ETAS-0 yields the best performing forecast, closely followed by ETAS-v (Figure 8 middle right panel).It not surprising that the versions of ETAS-f under-perform in this subregion as there should be no influence of the fluids.This also further supports the idea that the bad overall performance of ETAS-fv is due to the model interpreting the influence of fluid injection and/or production as variations in the background.
Hverahlíð sees a similar pattern, with ETAS-0 remaining the most well-suited model (Figure 7 bottom right panel).The area is relatively seismically active with 224 events recorded during the learning period, and 282 events are recorded during the forecasting period.The poorer performance of ETAS-f and ETAS-fv models in the subregion suggests that for production areas, the cumulative extracted volume is not a good proxy for the seismicity rate.Induced earthquakes in production areas (for oil and gas or geothermal energy) are often a response to reservoir compaction caused by long-term production.Hverahlíð has been in full operation since in late 2017, so these long-term effects are unlikely to be already visible.The area was also not associated with particular vertical displacement in any geodetic data and InSAR between 2015 and 2018 by Ducrocq et al. (2021).In the future, we will investigate other metrics to account for fluid production, for example, including a proxy term for poroelasticity and/or pore-pressure instead of a coupling with produced volumes per day.
In the reinjection areas, the results are more contrasted.Gráuhnúkar only sees moderate seismicity at rather smooth rates (Figure 6 bottom left panel); 99 events were recorded during the learning period, and 117 events are recorded during the forecasting period.The evolution of the relative performance of the Seismogenic Index is quite interesting, as it outperforms all models for the first half of the forecasting period, before taking a dive in July 2020 when an increase in the seismicity rate in the region happens (Figure 8 bottom left panel).In the Gráuhnúkar area, ETAS-0 and ETAS-v perform better than the versions with fluid forcing despite the injections being significant with roughly 12,000 m 3 /day injected (Figure 4b).The seismicity is known to be associated with the injections as the area was quiet beforehand (Flóvenz et al., 2015;Ritz et al., 2021).We would expect that ETAS-f picks up the injection-induced characteristics, however, the injections have been sustained since 2006 at very similar rates, which could explain how ETAS interprets the low-level seismicity as background, as our learning period doesn't cover pre-injection periods.
Húsmúli on the other hand, is the main seat of seismicity in the whole region; 579 events were recorded during the learning period, and 1,110 events are recorded during the forecasting period (Figure 6 middle left panel).
The reinjection in Húsmúli is mixed with production at the eastern border of the area, with compound volumes varying between 30,000 and 40,000 m 3 /day (Figure 4b).The area is known for its reactivity to injection changes (Gunnarsson et al., 2015;Ritz et al., 2021), which is confirmed in the models with ETAS-f outperforming all other models.The area also sees drastically varying rates of seismicity with three high-intensity events in October 2020, November 2020, and January 2021, when seismicity rates surpass 75 events per week.In this later part of the forecasting period, all advanced ETAS models see a stark change in their predictive capabilities and start performing better than ETAS-0 (Figure 8 middle left panel).These periods of high seismic intensity are the ones that see the Seismogenic Index's performance plummet (Figure 8 is cut-off at 600, but by the end of the forecasting period, the Seismogenic Index's information gain accumulates down to 3,600 in the Húsmúli subregion).
The top panel of Figure 8 shows the Seismogenic Index as the best performing model for the first half of the forecasting period.However, the only areas where this is highlighted in the subregions panels is Gráuhnúkar, where the performance is marginally better than all the ETAS models.This means that outside of our areas of interest (and in the cells where there is injections and/or production), the Seismogenic Index performs quite well compared to the ETAS models.

Adaptability to Rapid Seismicity Rate Changes
We now focus on the Húsmúli reinjection area that sees semi-periodic bursts of seismic activity which do not seem to be linked to changes in the injection rate, temperature or well-split (Figure 9a).The origin of these highintensity periods of seismicity is unclear, they could be induced, linked to the volcanic activity of the region or to the rift.In this section, we investigate how fast models are able to adapt to the new baseline seismicity during these high-intensity events.

ETAS Parameters Interpretation
The branching ratio η represents the share of events categorized as triggered aftershocks.Shifting from ETAS-0 to a varying background rate (ETAS-v) leads to more events being classified as background, thus the lower branching ratio (Figure 9c; Mizrahi et al. (2023)).Similarly, the ETAS-f model allows for more events to be classified into their two background categories (natural and fluid-induced), leading to an even lower base value of the branching ratio η.The effects of the free background and the additional fluid-background term add up to make ETAS-fv the model with the lowest branching ratio of the class.
The tectonic background term μ tect (given as number of events per km 2 per day in Figure 9) also shows this effect of increasing classification into 'background' with model complexity, with ETAS-0 at the lowest global values, ETAS-f with a significant increase in its base value due to the added fluid-background term which gives the model more flexibility to interpret events as background earthquakes, and the ETAS models with varying background (ETAS-v and ETAS-fv) at the highest base values of μ tect (mean value on the active area, Figure 9b) because they give the most flexibility to this term.If we look at the evolution of the background term in time, we notice a general decrease as the learning period length increases.Mizrahi et al. (2021b) proposed that this downwards trend denotes that longer learning periods allow models to reveal more long-term earthquake interaction and aftershock sequences, leading to more events being classified as "triggered" as the forecasting experiments progress.This general trend of μ tect is interrupted by bursts of seismic activity in the field, as highlighted by the gray areas marking the "high-intensity periods" in Húsmúli.All ETAS models see an increase in the background rate to accommodate this new temporary baseline seismicity rate.The relaxation of μ tect after the high-intensity periods is slow and seems to follow the general decreasing trend that we interpret as being due to the lengthening of the learning period.
During the high-intensity periods, the branching ratio for all models increases to accompany the increase in background seismicity rate.However, this increase is only temporary for the ETAS versions without fluid-forcing and η quickly recovers to its pre-high-intensity-period levels, whereas the ETAS-f and ETAS-fv models see a much slower recovery of η with a permanent change in the parameter resulting from the high-intensity period.These shifts in μ tect and η suggest that all ETAS models accommodate the high-intensity periods by producing more background events that are also produce more aftershocks.
For the ETAS with fluid forcing class, the fluid-background term μ f luid (also given as number of events per km 2 per day in Figure 9) combines the scaling factor ι (Equation 12) with cumulative compound volume.While ι defines the seismicity induced by a unit volume, μ f luid gives the rate of fluid induced events given the compound volume in the forecasting period.ETAS-fv shows higher base values of the fluid-background term, mirroring the higher values of μ tect .The evolution of μ f luid shows a non-systematic behavior (Figure 9d).The first and third high-intensity periods show an increase in the productivity of fluid driven events, which is what we expect ETAS-f and ETAS-fv models to do in order to adapt to increasing seismicity rates in an area with active injection and/or production.However, the second high-intensity period shows a sharp drop of μ f luid , maybe suggesting that ETASf models try to explain this change in the seismicity with tectonic events.When looking at the high-intensity periods on the map (Figure 9f), the first and third events cluster to the north where less fluids are distributed while the second high-intensity period clusters on the edge of the most heavy reinjection region, closer to the injection wells.This might suggest that peripheral regions can be more sensitive to fluid during certain periods, while the center of the area has so much fluid injected/extracted continuously that ETAS-f models can't explain a sudden increase of seismicity with a change of productivity of a unit of fluid, thus a drop in μ f luid down to its base value at the start of the second high-intensity period.
In the Ölkelduháls area the split background terms and set-up with the entire region used for training lead to the background terms in the case of ETAS-f being larger than for the standard ETAS-0.This calibration on the entire region also leads the triggering kernels to be different for ETAS-0 and ETAS-f, which drives the models to assign events as "background," "triggered" or "induced" differently.Thus, even in areas without fluids, ETAS-f provides a different forecast than ETAS-0.If the ETAS-f model had been trained only in Ölkelduháls, ETAS-0 and ETAS-f would yield the exact same forecast.Figure S2 in Supporting Information S1 shows the evolution of the parameters describing the aftershock distribution in time and space for all models of the ETAS class.

Seismogenic Index Parameter Evolution
The average value of the productivity parameter a f b increases in time to accommodate the change in seismicity rate (Figure 9e). Figure 10 shows a break down of the change of a f b in time in the three subregions where fluids are injected and/or produced.The activation feedback is much higher in the reinjection regions and Hverahlíð production area than the average shown in Figure 9e, highlighting the regions as seats of seismicity.The value of a f b is quite stable in Hverahlíð, where the production rate and seismicity rate are relatively constant during the forecasting experiment (Figures 4 and 7).On the contrary, in the injection regions of Gráuhnúkar and Húsmúli, where both the seismicity and injection rates have a more jagged behavior, the value of a f b sees sharp increases to try to accommodate the changing rate of seismicity.The case of Gráuhnúkar in particular shows that these "spikes" in a f b are followed by a slow healing phase during which the parameter is gradually lowered.
If we compare the fluid scaling factor ι for ETAS-f models and the Seismogenic Index's a f b as analogous "seismic productivity per unit volume" terms, we notice a similar relative increase in the first half of the forecasting period.However, the fluid productivity terms show different reactions to the high-intensity periods.The Seismogenic Index model only has this parameter to adapt to the fast changes in the seismicity rate, while ETAS-f and ETAS-fv models are able to explain the shift in seismicity rate with other parameters for example, assigned them as background or naturally triggered events.

Discussion
All the models used seem well suited to modeling seismicity in the Hengill geothermal field, although they show their respective limitations in specific subregions or following unexpected and abrupt changes in the seismicity rate.
However, the set-up of our pseudo-forecasting experiment also reveals some of its limits.First of all, the modeled period needs to match the length of operations in the field to provide a complete picture of the processes.As we saw in the Gráuhnúkar area, a 14 months training period on a reinjection in operation since 2006 with stable seismicity and injection rates can lead to a misattribution of the induced seismicity as background by ETAS-f and ETAS-fv models.The choice of the learning and forecasting periods was dictated by the high quality seismic catalog, however, it would be interesting to run a much longer term experiment that encompasses pre-reinjection times to see if ETAS-f and ETAS-fv models pick up on the onset of the injection-induced seismicity.
Following up on this point, the discrimination between induced and natural seismicity that ETAS-f models allow has not been fully explored.In particular, one could look at the attribution of aftershock sequences to induced or natural events.So far, the use of both ETAS and Seismogenic Index models only tell us broad information, for example, the case of Ölkelduháls, which must be a natural sequence (or induced by naturally occurring fluids) as no injection or production is taking place in the vicinity.
Another important aspect of induced seismicity mitigation methods relies on linking modeling results to hazard and risk calculations (Broccardo et al., 2019;Schultz et al., 2021Schultz et al., , 2022)).In this work, we stopped on the level of forecasting seismicity rates and did not venture further, however, we do plan to expand to hazard and risk in future work.

Seismogenic Index Model
In this study, we used a simple implementation of the Seismogenic Index, and could likely improve its performance in two ways.First, by improving the way how to account for the uncertainties of the model parameters replacing the Poissonian assumption by a method sampling from the parameter space around the optimized model parameters and defining a confidence interval around them (Rinaldi & Passarelli, 2021).Furthermore, as suggested by Broccardo et al. (2017), a Bayesian optimization approach using prior-distributions of the model parameters gathered from past induced seismicity cases in geothermal exploration sites, would be advantageous for two reasons: First, a model using the Bayesian theorem would automatically give a pseudo-Poissonian probabilistic output in form of a posterior distribution in terms of for example, seismicity rate.Second, such a model would be able to give a forecast also for cells that do not present an event during the training phase, contrary to the currently implemented version.A further way of improving the performance of the Seismogenic Index model is to improve the distribution of volumes in space, to take into account the field anisotropy driven by faults and fractures, as well as well directionality.Figure S3 in Supporting Information S1 shows an example of circular versus elliptical Gaussian distribution to account for well direction and faults in the Húsmúli region (Bessason et al., 2012;Khodayar et al., 2015).These results show that the distribution of volumes improves the performance of the Seismogenic Index by about 5% by accounting for a first order of subsurface anisotropy.

ETAS Models
In this work, we tested different variations of ETAS, from the standard implementation to a version with fluidforcing.However, this is only the beginning as many more things could be implemented into ETAS to increase its predictive skill in complex systems.For example, we saw that the ETAS with fluid forcing class seems promising for reinjection regions but performs poorly in production regions.Production induced seismicity is known to be linked to poroelastic stress changes (Segall, 1989;Zbinden et al., 2017).We want to introduce a proxy for poroelasticity in ETAS to account better for this phenomenon.This ETAS-poroelasticity would need to be tested separately on production regions, to emancipate ourselves from the noise of reinjection-induced seismicity and natural seismicity present in the Hengill field (Goebel et al., 2017;Segall & Lu, 2015).
In our current implementation, we distributed the volumes without accounting for the faults and geological unit adding anisotropy to the underground.We however have hydro-geological models at our disposal both for the Húsmúli reinjection area and Hengill-wide field (TOUGH2 models provided by Reykjavik Energy; Ritz et al. (2021)).These models could be coupled to ETAS (and to the Seismogenic Index model) into a hybrid model to give ETAS and the Seismogenic Index model a more realistic representation of the flow rate and volume distribution in the field.Such models with high-resolution local information on the fluid flow could help support operators in the well-sitting process.Indeed, the combination of a calibrated fluid-and heat-flow model with seismicity forecasting models provide a rare insight to decide on the development of production and injection areas while combining productivity information and a proxy for the seismic risk.
Besides the implementation of fluid-forcing, ETAS models can be made more flexible by playing on a number of parameters.In particular, adaptive smoothing kernels can be implemented in which the standard deviation of the kernel itself is determined to best fit the data (Nandan, Ram, et al., 2021;Ogata, 2011;Zhuang et al., 2002).Such models have been shown to have powerful forecasting capabilities.Furthermore, the productivity of aftershock sequences have been shown to have statistical differences between tectonic and induced earthquake sequences (Karimi & Davidsen, 2021;Llenos & Michael, 2013).One could account for these different signatures by having different parameters (e.g., productivity) for fluid-induced and tectonic events.

Other Models
Beyond ETAS and Seismogenic Index models, we could test other statistical models in a similar pseudoforecasting framework.For example, we could design an extremely simple model like "the rate that one will observe in the forecast is the same as the average in the last X weeks," and use it as a reference model instead of comparing all models to ETAS-0.
Purely statistical description of earthquakes in time and space are however limited as they do not account for the physical mechanisms behind induced seismicity.These mechanisms include static and dynamic stress transfer, the frictional properties of faults, as well as changes in pore pressure which are fundamental to explain induced events.Physics-based models able to describe aftershock triggering (e.g., Coulomb models; King et al. (1994); Catalli et al. (2016)) could be useful in complex contexts like Hengill.Physics-based model range from the simpler with analytical solutions for pore-pressure to more advanced ones like rate-and-state or models incorporating elements of poroelasticity.This last point would be interesting to compare in the production area of Hverahlíð where we saw that using the produced volume as a proxy is not very successful.Hardebeck (2021) notes a general under-performance of Coulomb-Rate-and-State models relative to ETAS models, but argues that incorporating heterogeneities in background conditions into physical forecasting models may be key in improving their performance.This conclusion is supported by Mancini et al. (2019), which however argues that stress-based models which consider secondary triggering mechanisms (stress changes, earthquake interactions) can preform similarly to ETAS in complex sequences.

Conclusion
The complexity of the seismicity in the Hengill geothermal field, with natural and induced events, as well as sudden periods of high-intensity activity, makes it an interesting but challenging testing ground for forecasting models.As we've seen, the Seismogenic Index and ETAS models have captured some seismicity aspects when forecasting the seismicity in the field and its different subregions, although the former can only produce a forecast in active injection and/or production regions.In such a complex geothermal field, the Seismogenic Index model performs similarly or better than ETAS models as long as the rate of seismicity remains relatively stable.As soon as an abrupt change in the seismicity rate occurs, the Seismogenic Index model perform significantly worse than any ETAS model.One major advantage of ETAS models lay in their ability to forecast natural seismicity, which allows them to fit data even in areas where no fluids are injected (as seen in Ölkelduháls).
Although there are still routes to explore to improve individual models, as highlighted in the Discussion section, the respective strengths of the models could be harvested by using an ensemble modeling framework.The concept to combine different models to obtain more robust forecasts and limit the individual model's biases, has been widely tested in recent years with all sorts of forecasting models (Bayona et al., 2021;Dempsey & Suckale, 2017;Király-Proag et al., 2018;Llenos & Michael, 2019;Marzocchi et al., 2012;Mizrahi et al., 2023).Such an approach with ETAS and Seismogenic Index-type models could prove useful as near-real-time tools in areas experiencing induced seismicity like Oklahoma or near the Groningen gas field in the Netherlands.
Hydraulic data is often difficult to stream in real time and might not be available with fine enough granularity to be useful for forecasting exercises.ETAS-0 does not require this input of fluids to yield forecasts.In the event that injection and/or production rates are available in real-time, ETAS with fluid forcing performs slightly better than standard ETAS and could easily be run in parallel to a Seismogenic Index model in an operational forecasting tool.

Figure 1 .
Figure 1.Location of the Hengill geothermal field in southwest Iceland (lower left panel) and distribution of seismicity in the Hengill geothermal field from 01.12.2018 to 31.01.2021,COSEISMIQ high quality hypo-DD relocated catalog (Grigoli et al., 2022).Dot size is proportional to the magnitude, color coded by time of occurrence.The dashed rectangles in the main panel and lower left panel show the active Hellisheiði geothermal area in which the models are evaluated.

Figure 2 .
Figure 2. Distribution of injection and production wellheads in the active geothermal area, as well as recorded seismicity between December 2018 and January 2021.The regions highlighted by boxes are areas of specific seismological interest.The grid lines show the space discretization for the forecasts.The direction of S Hmax indicated in the top left corner is based on borehole measurements in Húsmúli(Batir et al., 2012).

Figure 3 .
Figure 3. Statistical analysis of the catalog.(a) Joint estimation of the magnitude of completeness (vertical line) and b-value (blue errorbars) for different magnitudes; K-S distance (purple) and p Mc (yellow).(b) Frequency magnitude distribution of the catalog.The blue dots show the non-cumulative FMD, the vertical line denotes the magnitude of completeness, and the pink lines show the estimated b-value and its 1σ error.

Figure 4 .
Figure 4. (a) Distribution of cumulative compound volumes at the end of the experiment.(b) Temporal evolution of the compound volume in the areas of interest.

Figure 5 .
Figure 5. Schematic illustration of the pseudo-prospective forecasting experiments.

Figure 6 .
Figure 6.Comparison of the number of events recorded (dashed line) and forecast by the different models (colored lines).The top panel shows the entire Hengill field, the subpanels show the subregions.

Figure 7 .
Figure 7. Cumulative spatial joint log-likelihood in the active region (top panel) and subregions.

Figure 8 .
Figure8.Cumulative information gain of ETAS and Seismogenic Index models calculated relative to ETAS-0 in the active region (top panel) and subregions.The vertical scale is cut-off at 600, but the Seismogenic Index plunges all the way down to 2900 by the end of the forecasting period.

Figure 9 .
Figure 9. (a) Weekly number of events recorded in Húsmúli (moving average with 1 day increment; red) and injection rate in the area (blue).(b-d) Evolution of the ETAS models parameters.(e) Seismogenic index parameter evolution.The gray areas highlight the high-intensity periods.(f) Spatial distribution of the seismicity during the high-intensity periods overlaid on the cumulative volume map.

Figure 10 .
Figure 10.Evolution of the activation feedback in time during the forecasting experiment by subregion.