Introduction
Induced seismicity is associated with different energy technologies causing stress changes in the subsurface (for an overview see Foulger et al., Reference Foulger, Wilson, Gluyas, Julian and Davies2018 and National Research Council, 2013). Most if not all induced earthquakes occur on faults which are tectonically pre-stressed (Foulger et al., Reference Foulger, Wilson, Gluyas, Julian and Davies2018). In the 1960s, it was first recognised that earthquakes can be induced by injecting fluid into the subsurface (Healy et al., Reference Healy, Rubey, Griggs and Raleigh1968). These observations were explained by a model in which the effective normal stress on pre-existing fractures is reduced. By raising the fluid pressure, fractures are destabilised when the ratio between shear to effective normal stress exceeds their frictional resistance (Hsieh & Bredehoeft, Reference Hsieh and Bredehoeft1981). This is the most common mechanism to explain seismicity caused by fluid injection in geothermal reservoirs (Evans et al., Reference Evans, Zappone, Kraft, Deichmann and Moia2012; Ge & Saar, Reference Ge and Saar2022; Zang et al., Reference Zang, Oye, Jousset, Deichmann, Gritto, McGarr, Majer and Bruhn2014), by wastewater disposal (e.g. Schoenball et al., Reference Schoenball, Walsh, Weingarten and Ellsworth2018) and by hydraulic fracturing (Schultz et al., Reference Schultz, Skoumal, Brudzinski, Eaton, Baptie and Ellsworth2020).
Additional mechanisms were discussed that could be relevant for the occurrence of earthquakes in geothermal reservoirs. These include poro- and thermo-elastic stresses associated with fluid injection (Buijze et al., Reference Buijze, Bogert, Wassing and Orlic2019; Jeanne et al., Reference Jeanne, Rutqvist and Dobson2017; Rozhko, Reference Rozhko2010). The relevance of these mechanisms, however, remains unclear as it is difficult to relate observed earthquakes to a specific mechanism. In several geothermal reservoirs, fluid injection led to earthquakes causing material damage to buildings (Häring et al., Reference Häring, Schanz, Ladner and Dyer2008; Kim et al., Reference Kim, Ree, Kim, Kim, Kang and Seo2018; Schmittbuhl et al., Reference Schmittbuhl, Lambotte, Lengliné, Grunberg, Jund, Vergne, Cornet, Doubre and Masson2021). This repeated occurrence of damaging earthquakes brought the aspect of an induced seismicity risk into public and scientific focus.
In Europe, geothermal energy is regarded to play an important role for the sustainable energy mix in the future. A prerequisite for a successful growth of the geothermal industry is a safe extraction of geothermal energy. In our case example here, located in the Netherlands, this is explicitly stipulated by the Dutch State Mining regulator SodM (Jharap et al., Reference Jharap, van Leeuwen, Mout, van der Zee, Roos and Muntendam-Bos2020). In this context, the risk of induced earthquakes is one of the key factors that needs to be successfully managed. Several screening tools have been proposed to assess the potential that geothermal activities (or fluid injection in general) may induce seismicity (e.g. Davis & Frohlich, Reference Davis and Frohlich1993; Majer et al., Reference Majer, Nelson, Robertson-Tait, Savy and Wong2012) or that induced seismicity risk may become of concern for a specific project (Trutnevyte & Wiemer, Reference Trutnevyte and Wiemer2017). The current practise of seismic hazard assessment of geothermal projects in the Netherlands includes a Quick-Scan, which is a first and simple step to assess the project-specific potential for inducing seismicity based on a set of key operational and geological parameters (Baisch et al., Reference Baisch, Koch, Stang, Pittens, Drijver and Buik2016). Depending on the outcome of the Quick-Scan, a more detailed seismic hazard analysis (SHA) may be required.
As the shift towards new and environmentally clean energy systems, like geothermal plants, is continuously progressing in Europe, induced seismicity must be appropriately addressed in terms of hazard management. Hazard management procedures have been implemented and first executed by several countries including, for example, the United States (Wong et al., Reference Wong, Pezzopane, Dober and Terra2010) and Europe (Baisch et al., Reference Baisch, Carbon, Dannwolf, Delacou, Devaux, Dunand, Jung, Koller, Martin, Sartori, Secanell and Vörös2009). However, these seismic hazard assessment studies are rarely publicly available, not necessarily continuously updated and have not been methodically and scientifically evaluated. Case examples of repeated SHAs can provide a basis for such a methodological evaluation and reveal new findings about seismicity-inducing mechanisms and geological target formations.
The Californië geothermal field near Venlo is one of the rare examples providing a total of four SHAs conducted over a period of four years. This geothermal field near Venlo has provided heat for horticultural greenhouses for several years until a felt ML1.7 earthquake ultimately led to the suspension of geothermal activities in 2018. Here, we present the four seismic hazard assessment studies conducted for the Californië project, as well as an earthquake interpretation study and a hypocentre relocation study. We discuss how findings of these studies resulted in risk mitigation procedures and how successful and appropriate these risk mitigation procedures were in retrospective. In addition to that, we reveal that the newly inferred mechanism of thermal contraction provides a possible explanation for the induced seismicity observed at Californië.
The Californië geothermal project is the first project in the Netherlands to target the Dinantian, which is the main geothermal play considered in the national geothermal energy development strategy (Heijnen et al., Reference Heijnen, Jaarsma, Veldkamp and Schavemaker2019). Therefore, the induced seismic risk management at Californië is an important showcase for future geothermal project development in the Netherlands, but also in Germany, Belgium and possibly elsewhere.
Geothermal project Californië
The Californië geothermal system is located in the municipality of Horst an de Maas near Venlo in the south-east of the Netherlands, close to the German border (Fig. 1, left). Tectonically, the site is situated in the Venlo Block, a stable fault block, north-east of the Peel Block and Roer Valley Graben, which are part of the active Roer Valley Rift system. The Venlo Block is bounded by the Tegelen Fault in the south-west and the Viersen Fault Zone to the north-east (Houtgast & van Balen, Reference Houtgast and van Balen2000).
Fig. 1. Left: Location of the Californië geothermal project (red square) and main structural elements of the Roer Valley Rift System. Earthquake epicentres (KNMI-catalogue as of 2015 when the first SHA was conducted) are plotted as open circles scaled by magnitude (magnitude range M = −0.1 to M = 5.8). Right: Zoomed in section around the red square showing well trajectories of the CWG / CLG doublets (coloured lines) and fault trajectories at the top of the Carboniferous limestone group (black lines). Fault trajectories were derived from seismic interpretations and implemented into the numerical model accordingly. The major Tegelen fault extends farther than the seismic mapping indicates (dashed lines). Labelled lines indicate top reservoir depth level (in m). The grey-shaded zone denotes the extension of the Tegelen fault with depth. The Black arrow points in the Northern direction. RD (Rijks-Driehoek) coordinate system used. The red dashed line denotes the location of the vertical section as shown in Fig. 2. Modified figures from Vörös et al. (Reference Vörös, Baisch and Stang2015b).
The geothermal system consists of two doublets (Fig. 1, right). The first doublet was operated by Califonie Wijnen Geothermie b. v. (CWG) and started operations at the end of 2013. The second doublet (operated by Californië Lipzig Gielen Geothermie b. v., CLG) became operative in June 2017 and is located approximately 2 km to the North of the CWG doublet. The production wells target the Dinantian Carboniferous limestones at a depth below 2 km. Both wells produced hot fluid of around 75°C from the Tegelen fault for heating of horticultural greenhouses. The cooled down fluid was returned into the same formation through the reinjection wells further away from the fault (Fig. 2).
Fig. 2. Production (red, either CWG or CLG) and injection wells (blue: CWG, light blue: CLG) in a schematic geological cross-section from West to East through the project area. The location of the cross-section is outlined in Fig. 1 (red dashed line, not the same scale). The doublets are offset perpendicular to the drawing plane. The trajectories of the production wells penetrate the Tegelen fault. The injection well of the CWG doublet was originally drilled into the Tegelen fault but subsequently got blocked. Fluid was reinjected through a slotted liner into the reservoir formation at a depth range between 1800 and 2050 m (blue). The trajectory of the CLG reinjection well is directed away from the Tegelen fault. The geothermal aquifer comprises the Graben-related Tegelen fault zone and the Dinantian Zeeland formation (Kolenkalk Group). Modified figure from Vörös et al. (Reference Vörös, Baisch and Stang2015b).
The CWG doublet consists of the wells GT01 and GT03 and a third well (GT02). GT02 was initially planned as a reinjection well but abandoned due to a blockage. The production well GT01 intersects the Tegelen fault at 2100 m depth. In contrast to that, the reinjection well GT03 penetrates the reservoir layer at a shallower depth of approx. 1500 m and a distance of 450 m away from the fault. The GT03 well was originally drilled into the Tegelen fault and subsequently got blocked. Based on well logging tests, it is assumed that the lower well path is hydraulically isolated and that fluid injected into GT03 solely entered the Dinantian formation through a slotted liner at a depth between approx. 1380 and 1590 m.
The second geothermal doublet, CLG, consists of the production well GT04 and the injection well GT05. The production well intersects the Tegelen fault at 2390 m depth. The fluid was reinjected in GT05 at 1640 m depth into the reservoir. Reinjection occurred at a lateral distance of 2 km from the production well and the Tegelen fault.
Fluid was produced at rates of 400 m3/hrs (CWG) and 200 m3/hrs (CLG) with reinjection pressures of 1.2 and 2.2 MPa, respectively (Fig. 3).
Fig. 3. Measured flow rate (top) and wellhead pressure (bottom) at the reinjection wells GT03 (CWG, black) and GT05 (CLG, orange). The CWG doublet started production already at the end of 2013. Data, however, are only available after January 2014.The occurrence of induced earthquakes is denoted by black circles scaled by magnitude (top). The red dashed line denotes the occurrence time of the felt earthquake on 3 September 2018. The felt event is the largest event (large circle) and is the third event within a cluster of events between 3 September and 9 September. The timeline of seismic monitoring and studies performed related to the seismic hazard at the site is indicated at the top of the figure. Modified figure from Baisch & Vörös (Reference Baisch and Vörös2019).
Seismic monitoring of the CWG doublet commenced with three monitoring stations in September 2014 (stations K01, K02, and K03 in Fig. 7). Two additional stations (K04 and K05) were deployed in November 2015 before drilling of the wells for the second geothermal doublet.
Each monitoring station is equipped with a 3-C short-period surface seismometer (Lennartz LE3D) sampled at 100 Hz. Data are transmitted by cell phone modems to a data centre, where automatised processing is performed in real-time.
Since the start of seismic monitoring in September 2014, a total number of 17 local earthquakes in the magnitude range ML = −1.2 to ML = 1.7 have been detected. The earthquakes occurred between Aug 2015 and Sep 2018 in time intervals where one, both or none of the two doublets were in operation (Fig. 3). As will be discussed in more detail below, most of the earthquakes occurred after geothermal production had stopped or after the production rate was reduced. Although epicentre locations of the earthquakes correlate with the location of the CWG wells, initial estimates of the earthquake depth indicated a deeper origin between 4.5 and 6 km depth. These depth estimates were based on the range of seismic velocities typically assumed for the Netherlands. Since the region exhibits natural seismicity, the cause of the earthquakes could not immediately be pinpointed. Only at a later stage, a calibrated seismic velocity model could be used for re-locating the earthquakes at reservoir level (see section on hypocentre relocation).
Operations at the CWG doublet stopped in May 2018, as no permit for continuous operation was granted by the regulator. Production from the CLG doublet was suspended in August 2018 following an ML = 0.0 earthquake. Six days after stopping operations, a felt earthquake of magnitude ML = 1.7 occurred within the previous cluster of reservoir seismicity. This was followed by eight additional earthquakes ML ≤ 0.0 within six days. No further earthquakes have been detected until today.
Following the earthquake sequence, no operational permission was granted by the regulator. The future of the project is currently unclear.
Seismic hazard assessment
Induced seismic hazard assessment for the Californië geothermal system was a continuous process that started at an early stage of the project. Reviews and updates of the seismic hazard assessment were conducted to account for operational changes and the occurrence of seismicity.
Here, we summarise the main findings of the seismic hazard assessments as originally published. We maintain the level of information available at the time when the hazard assessments were conducted, using original diagrams with editorial modifications only.
Keeping the original content and preserving the chronological order provides the basis for evaluating the induced seismic risk management at Californië. The evaluation follows in the discussion section and is based on all available information.
Although the Dutch ‘framework for seismic hazard assessment in geothermal projects’ (Baisch et al., Reference Baisch, Koch, Stang, Pittens, Drijver and Buik2016) was implemented only after the first geothermal doublet at Californië had started production, the induced seismic risk management at Californië was largely consistent with these guidelines. An initial Quick-Scan (Baisch et al., Reference Baisch, Koch, Stang, Pittens, Drijver and Buik2016) indicated a medium potential for inducing seismicity, and a more detailed hazard assessment was required. The seismic hazard assessment was conducted for each doublet separately and was updated in the course of the project, as the exploitation scheme changed and more data on subsurface conditions became available. The following sections summarise the four seismic hazard assessments and two related studies.
SHA#1 (CWG, June 2015)
The seismic hazard assessment followed a physics-based approach, where stress changes associated with subsurface activities were numerically simulated as described in more detail below. Subsequently, a scenario technique was applied to study the seismicity response to simulated stress changes assuming two different scenarios of subsurface conditions. These scenarios reflect the uncertainty of subsurface conditions. An ‘expected scenario’ is defined based on the geothermal project developer’s interpretation of subsurface conditions. The second ‘extreme scenario’ considers unfavourable subsurface conditions, where the Tegelen fault behaves seismogenic and is close to stress criticality at reservoir depth.
Quantifying the likelihood for the two individual scenarios was considered impossible. Instead, risk mitigation measures in terms of a traffic light protocol were defined such that damage-relevant seismicity should be avoided even in the ‘extreme scenario’. This strategy differs fundamentally from the concept of a probabilistic seismic hazard assessment (Cornell, Reference Cornell1968; McGuire, Reference McGuire1995), where risk and uncertainties are estimated quantitatively. Given the lack of previous (induced) seismicity at Californië and accounting for the non-stationary nature of induced seismicity, a probabilistic quantification of seismic hazard was considered not feasible.
The workflow of the seismic hazard assessment can be summarised as follows:
In a first step, the geological-tectonical situation at the project site was investigated (Figs. 1 and 2). Using a local interpretation of fault trajectories determined from 2-D seismic lines and the regional stress field model of Hinzen (Reference Hinzen2003), slip tendencies ST, defined as the ratio of shear to effective normal stress (Moeck et al., Reference Moeck, Kwiatek and Zimmermann2009), were calculated on all patches of the mapped faults. The resulting slip tendencies fall in the range from 0.9 to 1.1. Slip tendencies are higher than typical values of the coefficient of friction, indicating that the mapped faults in the subsurface are likely to be critically stressed. For example, for the reservoir rocks we assumed a coefficient of friction of μ = 0.6, which is in the typical range for limestones (Beblo et al., Reference Beblo, Berktold, Bleil, Gebrande, Grauert, Haack, Haak, Kern, Miller, Petersen, Pohl, Rummel and Schopper1982, table 11).
For the ‘expected scenario’, a geometrically simple reservoir model was implemented as a numerical 3-D Finite Element model in Comsol Multiphysics for simulating hydraulic pressure changes associated with (the planned) geothermal production. A Darcy flow regime (Bear, Reference Bear1975) was assumed. Poro-elastic stress changes were considered to be of secondary order and were thus discarded. Mechanical stress changes associated with thermal contraction were simulated separately using the semi-analytical solutions by Okada (Reference Okada1992) following the approach described in Baisch et al. (Reference Baisch, Carbon, Dannwolf, Delacou, Devaux, Dunand, Jung, Koller, Martin, Sartori, Secanell and Vörös2009) and parameters as listed in Tables 1 and 2. This approach approximates thermo-elastic stresses on faults located in the far field, that is beyond the cooled zone around the injection well. We have compared modelling results from the semi-analytical approach to those of a thermo-hydraulic-mechanically coupled finite element model and found a first-order consistency (see appendix).
Table 1. Parameter for the hydraulic model used in the two SHA updates. The model includes results from drilling and testing the CLG wells, resulting in a more detailed resolution of the subsurface with modified hydraulic parameter (layers of the Zeeland group are labelled L2-L5). Layers are thinning out in the eastern direction, and the two values denote the respective maximum / minimum thickness.
Table 2. Parameter for the computation of thermal stresses.
Numerical simulations indicate that after 5 years of production, temperature changes >0.5°C are limited to distances <325 m around the injection well and do not reach the Tegelen fault. Thermo-elastic stresses, however, can still be relevant for the seismicity evolution of the Tegelen fault (e.g. Jeanne et al., Reference Jeanne, Rutqvist and Dobson2017).
Stress perturbations on faults are described by Coulomb stress changes ΔCS (Scholz, Reference Scholz2002):
$${\rm{\Delta }}CS = {\rm{\Delta }}\tau - \mu \cdot \left( {{\rm{\Delta }}{\sigma _n} - {\rm{\Delta }}{p_{fl}}} \right)$$
with
$\delta \tau$
,
$\delta \sigma _{n}$
and
$\delta p_{fl}$
denoting changes of shear stresses, normal stresses and fluid pressure, respectively. μ denotes the coefficient of friction. Positive values of ΔCS shift the state of stress on a fault (patch) closer to failure. Coulomb stress changes along the Tegelen fault resulting from thermal contraction after 2 years of circulation at 250 m3/hrs are shown in Fig. 4. For this simulation, a pure normal faulting regime was assumed. The stress regime, however, is not well constrained, and a strike-slip regime was considered equally likely in the seismic hazard assessment. At a later stage, we will show that a similar pattern of Coulomb stress changes rotated by 90° results when assuming a strike-slip regime (Fig 10).
Fig. 4. Simulated Coulomb stress changes for a normal faulting regime related to thermal contraction on the Tegelen fault after 2 years of continuous circulation of the CWG doublet. A circulation rate of 250 m3/hrs has been assumed. Stress magnitudes are colour-scaled. The well trajectories are depicted as red (GT01, production) and blue (GT03, reinjection) lines. In the SHA, normal faulting as well as strike-slip faulting have been considered, both resulting in comparable stress magnitudes but different spatial patterns. The black arrow denotes the Northern direction. Coordinates with respect to x = 204,042, y = 380,050 (RD). Modified figure from Vörös et al. (Reference Vörös, Baisch and Stang2015b).
For the ‘extreme scenario’, an ad hoc assumption was made that already smallest stress perturbations of
$\Delta CS \ge $
0.06 MPa may trigger seismicity. Furthermore, it was assumed that the fault is half a stress drop away from stress criticality. Based on the numerical simulations of stress changes (see above), the scalar seismic moment of induced seismicity was estimated from the continuous fault area subjected to stress load with
$\Delta CS \ge $
0.06 MPa. Seismic moment was converted to earthquake magnitude using the magnitude definition of Hanks & Kanamori (Reference Hanks and Kanamori1979).
Since no seismicity was detected during the previous 1.5 years of geothermal production, it was considered unlikely that fluid overpressure would cause seismicity in the future if the geothermal doublet was operated at the same production rate. This conclusion was based on the numerical model (see appendix) indicating that quasi-stationary hydraulic conditions were reached already after 2 weeks of continuous geothermal production. It should be noted that the numerical model did not account for the temperature dependency of the fluid viscosity. The impact of temperature-dependent fluid parameters, however, was assumed to be small compared to, for example, the observed injectivity increase of the injection well.
While fluid overpressure was not considered a relevant factor for causing seismicity, thermo-elastic stresses on the Tegelen fault were identified as a potential issue if the rocks at reservoir level deform seismically. These stresses could result from the cold water injection close to the Tegelen fault and accumulate over time. The simulated seismicity response for the ‘extreme scenario’ shows a systematic increase of the maximum magnitude over time, reaching the level of Mw = 2 after approximately 2 years of continuous operation (Fig. 5).
Fig. 5. Potential magnitude increases with time related to Coulomb stress changes on the Tegelen fault (compare Fig. 4). It is conservatively assumed that a continuous patch on the fault, subjected to stress perturbations above the critical threshold of ΔCS = 0.1 MPa, fails simultaneously. Figure from Vörös et al. (Reference Vörös, Baisch and Stang2015b).
Thermo-elastic stresses simulated for other mapped faults were much smaller due to their larger distance from the injection well. The level of
$\Delta CS \ge $
0.06 MPa is only reached on the Tegelen fault.
For limiting earthquake strength, real-time monitoring with a traffic light protocol (TLP, Bommer et al., Reference Bommer, Oates, Cepeda, Lindholm, Bird, Torres, Marroquín and Rivas2006) was recommended as a risk mitigation measure (Table 3). The TLP was designed to prevent the occurrence of an earthquake causing even slightest material damage. The metric of the proposed TLP was peak ground velocity (PGV), which was related to material damage through engineering standards (SBR, 2010). A ‘red light’ threshold value of PGV = 0.3 mm/s was thus defined for stopping operations. This threshold corresponds to the lower level of human perceptibility and is by a factor of 10 lower than the level at which damage to (sensitive) buildings is considered possible (SBR, 2010). Implementing such a large safety margin was deemed necessary to account for typical trailing effects, where the strongest earthquakes occur post-operation (see Baisch et al., Reference Baisch, Koch and Muntendam-Bos2019 for a more recent discussion).
Table 3. Traffic Light Protocol (TLP) for the Californië geothermal system. The TLP is based on peak ground velocities (PGV) measured at the surface. A stop of operations (‘red light’) is triggered in case ground vibrations exceed the level of human perceptibility. This threshold value accounts for a potential increase of earthquake strength after stopping operations (‘trailing effect’) such that minor damage to building is avoided.
The injectivity of the GT03 well had systematically improved during the 1.5 years of previous geothermal production, indicating that flow paths in the reservoir may have changed. For example, fracture permeability may have improved locally due to fracture cleaning. To ensure the applicability of the hazard assessment, a pressure monitoring focusing on increasing pressure at constant flow rates was suggested as an additional mitigation scheme. Furthermore, a re-evaluation of the seismic hazard was recommended in case of new data and/or if operational parameters deviate from the range considered in the seismic hazard assessment.
SHA#2 (CLG, August 2015)
After the first geothermal doublet (CWG doublet) went into operation, planning for a second geothermal doublet nearby was commenced by a different operating company (CLG doublet). Extending the geothermal system was considered a delicate subject as stress interferences of the two doublet systems were expected due to their proximity. Prior to drilling the CLG wells, a seismic hazard assessment was performed for the CLG doublet. This SHA also accounted for stress contributions from the existing CWG doublet (Vörös et al., Reference Vörös, Baisch and Stang2015a).
The seismic hazard assessment followed the same approach as outlined in the previous section. Different geological interpretations of the two operators were synchronised within a simplified numerical model for simulating hydraulic pressure changes associated with the simultaneous operations of the two geothermal doublets.
Fig. 6 shows hydraulic overpressure for a production scenario in which both doublets are operated simultaneously. Pressure changes with
$\Delta p_{fl}$
= 0.1 MPa (equivalent to the assumed level above which seismicity can be triggered, Equation 1) are not observed at the Tegelen fault, but at the Velden fault to the north-east. Given the low level of overpressure, the potential for causing seismicity by pressure changes on the Velden fault was still considered low.
Fig. 6. Numerical model (left) and simulated fluid pressure changes Δpfl in the reservoir layer for a production scenario with the two doublet systems (right). Grey surfaces (left) show the aquifer and the Tegelen fault of the numerical model. Pressure changes shown after 190 days of continuous operation at a rate of 310 m3/hrs. Isobars are depicted in red, reservoir faults in grey, main faults are labelled. Arrow indicates Northern direction. Coordinates with respect to x = 204,042, y = 380,050 (RD). Modified figure from Vörös et al. (Reference Vörös, Baisch and Stang2015a).
Simulation of the stress changes caused by thermal contraction resulting from CLG production yielded minor stress changes on the mapped faults. These changes were considered to be too small to cause seismicity.
Incorporation of the cooling stresses resulting from the CWG injector (compare previous section) and the possibility of an unmapped, critically stressed fault in the vicinity of the CLG injector, resulted in the recommendation of seismic real-time monitoring with a TLP as a risk mitigation measure. The proposed TLP was based on the same threshold values as for the CWG doublet system (Table 3), while a red stoplight now implied the shutdown of both doublets.
SHA#3 (CLG update, March 2018)
The seismic hazard assessment for the CLG doublet was updated after drilling and testing the CLG wells (Vörös et al., Reference Vörös, Baisch and Stang2015b). Subsurface stress perturbations were simulated utilising a revised reservoir model with a modified geometry. Hydraulic conductivities had to be reduced to stay consistent with well-testing results. Consequently, simulated overpressure on the Velden fault slightly increased in the updated simulations. The simulations revealed that quasi-stationary conditions had already been reached during the previous well tests, which were not associated with measurable seismicity. Therefore, overpressure causing seismicity was considered unlikely if circulation rates were not to be increased.
Thermo-elastic stresses on the mapped faults (located beyond the cooled zone) were numerically simulated using the approach outlined in SHA#1. Sensitivity tests have shown that thermal contraction stresses are relatively insensitive to the geometrical details of the cooled zone as long as the faults are sufficiently far away from the cooling front. Therefore, it was concluded that the modified reservoir model had no significant impact on the outcome of previous simulations of stress perturbations resulting from thermal contraction.
Earthquake interpretation study, February 2019
Following an earthquake sequence in September 2018, a detailed study was conducted to investigate the cause of the earthquakes and their relation to the two geothermal doublets (Baisch & Vörös, Reference Baisch and Vörös2019). Hypocenter locations were determined with a linearised inversion scheme based on P- and S-phase arrival times (Baisch et al., Reference Baisch, Bohnhoff, Ceranna, Yimin and Harjes2002). Epicentre locations determined with the initial (uncalibrated) seismic velocity model clustered near the CWG doublet with an isolated earthquake near the CLG injection well (Fig. 7). Hypocentral depth was estimated around 6 km. This depth was interpreted to result from the uncalibrated velocity model. Based on the immediate response of the seismicity to production rate changes and excluding natural causes, it was speculated that the earthquakes were more likely located shallower, in the depth range from 2.5 to 3.0 km and thus within the range of stress perturbations from the geothermal system.
Fig. 7. Absolute hypocentre location in map view. Error bars show the location accuracy with a 2σ confidence level (formal inversion error). Triangles denote the location of the seismic monitoring stations. Events were colour-coded according to the time of occurrence (see legend). The grey patch depicts the Tegelen fault. The magenta lines show the well trajectories of GT01, GT03, GT04 and GT05, respectively. Black arrow indicates Northern direction. Coordinates with respect to x = 204,042, y = 380,050 (RD). The first six events occurred prior to production start of the CLG doublet. Event hypocentres are depicted as determined for the SHA before the modification of the velocity model. Figure from Baisch & Vörös (Reference Baisch and Vörös2019).
To evaluate the stress impact of the two doublet systems separately, the six earthquakes occurring prior to CLG operation were investigated. It was noted that these earthquakes systematically occurred after the CWG production rate was either reduced or production was stopped. This hypothesis was tested using a simple mathematical model comparing the production rate Q(t i ) at the time t i when the i’th earthquake occurred with the average production rate QAV i prior to the earthquake. The mean flow rate prior to earthquake i was defined as
$${\rm QAV_i}\left( T \right) = \;{{\mathop \int \nolimits_{ti - T}^{ti} Q\left( t \right)dt} \over T}$$
with T denoting the time window length over which the flow rate is averaged. Furthermore, the parameter
$${k_i}\left( T \right) = {{Q\left( {{t_i}} \right)} \over {\rm QAV_i}\left( T \right)} $$
was defined to compare the production rate at the time of the earthquake to the previous, averaged production rate. To avoid bias from production rate variations close to zero, production data were rounded towards 10 m3/hrs bins. Moreover, the numerical resolution limit (10−16) was added to the average rate QAV i (T) to obtain k i = 0 in case Q i = 0 and QAV i (T) = 0. If earthquake timing truly correlates with negative production rate changes, k i is systematically smaller than 1. Indeed, this was found for the 6 earthquakes occurring at a time when only the CWG doublet system was operated. Fig. 8 shows k i as a function of window length T. Two events occurred during shut-in (i.e. Q i = 0). These events exhibit a k i = 0, independent of T. For the remaining four events, k i is close to 1 for small T (approximately < 1.2 days) reflecting the delay time between rate changes and earthquake occurrence. For larger T, k i is generally smaller than 1.
Fig. 8. Ratio between production rate Qi at the occurrence time of seismic event i and the average rate QAV prior to the event as a function of time length T over which the production rate is averaged. The ratio is shown for all six seismic events occurring between 31 August 2014 (begin of seismic monitoring) and 31 May 2017 (begin of CLG production). Earthquakes occurring during shut-in (Qi = 0) show up as a flat line. Note: In a diffusion-type triggering model, delay times can vary even if all earthquakes occurred at the same location. Delay times are sensitive to the specific triggering level of each event. Modified figure from Baisch & Vörös (Reference Baisch and Vörös2019).
In a subsequent step, synthetic earthquake catalogues were used to investigate the probability of a coincidental correlation between production rate changes and earthquake occurrence time. 100,000 synthetic catalogues were compiled, each containing six earthquakes with an occurrence time sampled from a uniform random distribution. The rate at which the synthetic catalogues exhibit k i (T) smaller than observed was found to be in the order of 1e−3 and 1e−4. Therefore, a coincidental correlation between earthquake occurrence and CWG activities can be excluded at a high confidence level. Consequentially, changes of the hydraulic fluid pressure Δp fl were identified as the triggering mechanism for the 6 seismic events under consideration.
A conceptual explanation was proposed. In this model, the Tegelen fault was loaded by Coulomb stresses resulting from thermal reservoir contraction due to cold water injection close to the fault. These stresses were interpreted to be the root cause for the earthquakes. During production, fluid pressure on the Tegelen fault was lowered. This resulted in higher effective normal stresses
$\Delta\sigma _{n,eff} = \Delta \sigma _{n}- \Delta p_{fl}$
, stabilising the fault. At the lowered pressure conditions during production, the Tegelen fault was further loaded by additional thermal contraction stresses without becoming overcritical. In the proposed model, stress criticality and seismic failure occurred post-production when the in situ pressure in the Tegelen fault returned to equilibrium state. At that time, effective normal stresses were decreasing (thus destabilising the fault), while the additional stress load from thermal contraction remained.
Simple numerical simulations were performed to test the validity of the proposed mechanism. These simulations focussed on the triggering process and the associated timing of the earthquakes in response to changes in production history rather than attempting to model absolute pressure levels.
For example, Fig. 9 shows a simulated pressure evolution inside the Tegelen fault at 5 km depth. The timing of the earthquakes coincides with the sharp pore pressure build-up following shut-in. The level of pressure changes is very small in this model. At reservoir depth, pressure changes are larger but still imply earthquake triggering at a very low level of stress changes.
Fig. 9. Evolution of simulated fluid pressure changes in the Tegelen fault at a depth of 5 km. This conceptual model aims to explain the mechanism causing induced seismicity during time periods of reduced rates or shut-in. The model consists of a simple reservoir layer intersecting the Tegelen fault. Occurrence times of the seismic events (18-Aug 2015, 5-Dec 2015 and 26-Jan 2016) are marked by black arrows. Figure from Baisch & Vörös (Reference Baisch and Vörös2019).
The proposed model provided an explanation for the earthquakes occurring prior to the operation of the second doublet. The simultaneous operation of both doublets after September 2017 impeded an establishment of causal relations for the later events. It was noted, however, that the subsequent earthquake (occurring 8-4-2018) exhibited the same correlation with rate reduction at CWG while both doublets were operational.
The sequence of earthquakes in September 2018 occurred after the shutdown of both geothermal doublets. Relative hypocentre locations were determined indicating that these earthquakes occurred almost at the same location where the previous seismicity related to CWG production had occurred. By combining the observed spatial and temporal correlations, the study concluded that all but one earthquake was likely related to production at the CWG well.
The remaining earthquake (M = 0, 25-8-2018) was located further North in the vicinity of the CLG injection well. It was concluded that this earthquake could have been caused by CLG production but contributions from CWG production were also considered possible. With hindsight, it should be noted that the earthquake in question was shifted onto the Tegelen fault after re-locating hypocentres (see section on Hypocenter Relocation).
SHA CLG update, March 2019
As part of the permit application process for resuming production of the CLG doublet, the seismic hazard assessment for the CLG doublet was updated (Vörös & Baisch, Reference Vörös and Baisch2019). The study focused on the seismic hazard associated with resuming production at the CLG doublet with a decommissioned CWG doublet. This update followed the same approach used in the previous studies accounting for the seismicity observations and interpretations outlined in the Earthquake interpretation study (previous section).
It was concluded that resuming geothermal production at the CLG doublet will most likely not cause seismicity on the known faults. This conclusion was drawn based on the assumption that injection into the GT03 well, identified as the root cause for seismicity, will not resume.
The same TLP used in previous studies (Table 3) was proposed for limiting the strength of earthquakes potentially occurring on an unmapped fault. Given that the cause for one of the earthquakes was considered unclear (see previous section), it was recommended to stop production if an earthquake, independent of its magnitude, is detected near the CWG injector. The reasoning was that such an earthquake would indicate that seismicity can be triggered by extremely small stress perturbations, smaller than considered in the study. In that case, a re-assessment prior to recommencing production was deemed necessary as such extremely small stress changes were not taken into account in the study.
With these mitigation measures in place, the mitigated risk was considered to be acceptable.
Hypocenter relocation, January 2021
In a recent study, we have re-located hypocentres based on a seismic wave velocity model calibrated with recordings of several perforation shots fired in the well GT04 (de Pater, Reference de Pater2021). This perforation shot data were not considered previously as the shot time was not measured. To account for an unknown shot time, we have designed a new calibration procedure based on differential travel times between S-wave (t S) and P-wave (t P), which are independent of the shot time. Let v P and v S denote the compressional and shear wave velocities, and Δ the source-receiver distance, then
$$t_S - t_P = t_0 + \Delta /v_S - t_0 - \Delta /v_P = \Delta /v_S - \Delta /v_P$$
Rearranging yields
$$\Delta = v_P*v_S/(v_P - v_S)*(t_S - t_P)$$
In our approach, the mixed velocity term v P·v S / (v P-v S) was calibrated and subsequently used for re-locating the earthquakes using a linearised inversion (e.g. Baisch et al., Reference Baisch, Bohnhoff, Ceranna, Yimin and Harjes2002).
Re-located hypocentres aligned along the Tegelen fault close to the CWG wells approximately at reservoir depth (Fig. 10). Interestingly, also the most northern earthquake was consequently shifted towards the Tegelen fault. The cause of this earthquake, however, remained speculative.
Fig. 10. Re-located seismicity, based on the reviewed velocity model (black dots) and re-simulated thermal contraction stress changes ΔCS on the Tegelen fault. The view from South-West (left) and from top (right) is shown. Stress changes are colour-scaled in MPa according to the colour bar. The accumulated stresses have been re-simulated for May 2018, after the CWG doublet stopped operations. Simulations were based on the actual flow rate. Stresses were computed assuming a strike-slip failure regime. The resulting stress pattern coincides spatially with the re-located hypocentres of the events. The colour map is saturated at a value of 0.1 MPa. Maximum stress values of >0.3 MPa are obtained locally. Trajectories of the GT01 and GT03 wells are depicted as red and blue lines, respectively. The Northern direction is indicated by a black arrow. Coordinates with respect to x = 204,042, y = 380,050 (RD).
Discussion
In this article, we have documented the seismic hazard assessment performed for the Californië geothermal field in the Netherlands. Consecutive observations of induced seismicity resulting from geothermal operations were presented. For the evaluation of the seismic hazard assessment and the underlying geomechanical concepts, the depth at which the earthquakes have occurred is a crucial piece of information. This becomes evident when considering the mismatch between initial earthquake locations and the depth at which the largest stress changes were predicted. In our showcase, earthquake depth could be constrained only at a very late stage of the project. Dedicated calibration shots for calibrating seismic wave velocities could greatly improve seismicity interpretations in future geothermal projects.
Evaluation of seismic hazard assessment
The induced seismic hazard assessments conducted for the Californië geothermal field reflect the state of knowledge about subsurface conditions at different times. In this respect, the hazard assessments must be considered living documents, which is a complicating factor when evaluating their accuracy. To resolve this, we adopt the perspective of the regulator asking whether the hazard and risk level was sufficiently assessed at all times of the project and whether sufficient risk mitigation measures were always implemented to ensure safe operations. We base our evaluation on the induced seismicity observations.
With hindsight, it appears that the dominating processes controlling induced seismicity at the Californië geothermal field have been identified already in the first seismic hazard assessment. In this study, hydraulic and thermo-elastic stress changes were considered the most relevant factors driving the seismic hazard. While hydraulic pressure increase is the most common triggering mechanism considered for fluid injection scenarios, the (far field) stress impact resulting from thermal contraction has been discussed theoretically (e.g. Baisch et al., Reference Baisch, Carbon, Dannwolf, Delacou, Devaux, Dunand, Jung, Koller, Martin, Sartori, Secanell and Vörös2009; Jeanne et al., Reference Jeanne, Rutqvist and Dobson2017). However, case examples demonstrating the actual relevance of this process had been missing to this point.
Earthquakes observations at Californië are not consistent with the typical concept of post-injection pressure diffusion causing seismicity at greater distances from an injector (e.g. Hsieh & Bredehoeft, Reference Hsieh and Bredehoeft1981). At Californië, seismicity tends to occur within hours to a few days after a rate reduction or production stop (compare Figs. 3 and 8) and does not exhibit an apparent migration pattern with time (Fig. 7). The correlation between operational changes and earthquake timing indicates that the earthquake triggering process might be controlled by fluid pressure changes. This seems to contrast with the conclusion of the first seismic hazard assessment, where thermo-elastic stresses were considered as the key mechanism driving the seismic hazard. The subsequent earthquake interpretation study, however, provides an explanation for this correlation and at the same time supports the initial conclusion. At the Californië site, both doublets produce fluid from the Tegelen fault. During production, the highly permeable Tegelen fault is locally strengthened by fluid pressure reduction, serving as a ‘pressure shield’. Decreasing or stopping production removes this ‘pressure shield’ and subsequent pressure build-up can cause elevated fluid pressure levels locally, even with respect to undisturbed conditions. At the same time, thermal contraction stresses on the Tegelen fault accumulate, thus destabilising the fault. At those locations, where stress levels are close to failure, the increasing fluid pressure associated with the reduced production rates constitute the immediate cause of induced seismic failure. The root cause, though, is the steadily increasing thermal contraction stresses on the fault.
In this interpretation, we expect earthquakes to occur in those regions of a seismogenic fault, where thermo-elastic Coulomb stress changes are largest. In the initial seismic hazard assessment, these stress changes were simulated assuming a normal faulting stress regime (Fig. 4), while noting that a strike-slip regime is equally likely. With hindsight, we have simulated thermo-elastic stress changes for strike-slip failure on the Tegelen fault (Fig. 10). These simulations are based on the actual flow rate measured in the Californië wells. We note that most hypocentres are located close to the region of maximum Coulomb stress changes, thus adding further confidence in our interpretation.
We also note that the magnitude of the strongest ML1.7 earthquake, which occurred after injection of 8.9 million m3 of (cold) fluid into the GT03 well, is only slightly below the maximum magnitude of ML2.0 simulated in the initial seismic hazard assessment.
In this respect, the ‘extreme case scenario’ of the initial seismic hazard assessment has reasonably anticipated location and strength of the induced seismicity. This emphasises the value of a physics-based seismic hazard assessment for geothermal reservoirs, even when conducted with limited knowledge of subsurface conditions.
Inherent to the seismic hazard assessments conducted for Californië is the definition of a response protocol for avoiding a damaging earthquake even in the ‘extreme case scenario’. This is further discussed in the next section.
Earthquake timing and TLP
The design criterion for the TLP was to avoid material damage even when accounting for post-injection seismicity of comparatively large magnitude. In the hazard assessments, this ‘trailing effect’ was conceptually attributed to the process of post-injection pressure diffusion originating from the injection wells. The same concept has been proposed elsewhere to explain post-injection seismicity (e.g. Hsieh & Bredehoeft, Reference Hsieh and Bredehoeft1981, Baisch et al., Reference Baisch, Vörös, Rothert, Stang, Jung and Schellschmidt2010). At Californië, however, timing and location of the earthquakes are not consistent with this simple concept. This raises the question of whether the TLP at Californië still provided an efficient risk mitigation measure. We note that the TLP as designed in the seismic hazard assessment was not strictly tested since operations were suspended already 6 days prior to the occurrence of the stoplight earthquake.
In our interpretation, earthquake timing is dependent on the combined effect of thermo-elastic stress load and fault stabilisation due to pressure depletion. Cold water injection into GT03, which was interpreted as the cause for thermo-elastic stress load on the Tegelen fault, had stopped already approximately 4 months prior to the M1.7 earthquake. We interpret the delayed occurrence of the earthquake to result from the ongoing GT04 production, stabilising the Tegelen fault. The M1.7 earthquake occurred only after the GT04 well stopped production. In this scenario, earthquakes predominantly occur post-operational, thus questioning the concept of a TLP in general. A possible strategy for mitigating this type of trailing effect is regular interruptions of geothermal production to allow for pressure build-up and for testing fault stability. Such interruptions were not foreseen in the hazard assessments. We note, however, that the estimates of the (unmitigated) maximum earthquake magnitude performed in the hazard assessment studies are not depending on the performance of the TLP.
Conclusions
The Californië geothermal field provides a rare showcase for evaluating induced seismic risk management. Several seismic hazard assessment studies were performed for the geothermal doublets at Californië. These studies were based on deterministic models, where subsurface stress changes resulting from geothermal activities were numerically simulated. A scenario technique with an ‘extreme’ and an ‘expected’ scenario was used to account for uncertainties of subsurface conditions. Besides the ‘expected case’, the ‘extreme case’ of a critically stressed, seismogenic target fault was considered. A TLP was defined aiming to avoid damage-relevant seismicity even in the ‘extreme case’. Strength and location of subsequent earthquakes are reasonably consistent with the numerical simulations of the ‘extreme case’. In hindsight, we consider the risk level to be sufficiently assessed during all times of the project.
In contrast to similar fluid injection operations, our analysis indicates that fluid overpressure may not be the driving force for the seismicity at Californië. Instead, numerical modelling results suggest that thermo-elastic stresses resulting from cold water injection close to the Tegelen fault could be the root cause for the observed seismicity. These stresses accumulate over the lifetime of the geothermal system and can cause seismicity even beyond the cooled zone. The proposed mechanism could also explain the intriguing observation that most of the earthquakes occurred within hours to a few days after geothermal production had stopped or after the production rate was reduced. During production, thermo-elastic stress loading on the Tegelen fault is temporarily counterbalanced by pressure drawdown near the production wells. Following a rate reduction or a production stop, pressure build-up around the production wells destabilises the fault and might finally control the timing of earthquakes. For the Californië geothermal field, we speculate that re-locating the first injection well further away from the Tegelen fault would most likely eliminate the cause for the previous seismicity.
Acknowledgements
We would like to thank Pieter Wijnen (CWG) and Lodewijk Bourghout (CLG) for allowing us to publish the seismic hazard assessment and interpretation studies.
We gratefully acknowledge comments by three anonymous reviewers which helped to improve the initial manuscript.
Appendix
Simulation of hydraulic pressure changes
Different reservoir models were implemented within the COMSOL Finite Element software package (COMSOL Multiphysics®, www.comsol.com). For the initial SHA of the CWG doublet system, an existing reservoir model developed by the operator was numerically implemented. The model is based on 2-D seismic data interpretation, hydraulic testing and well logging data. It consists of the geothermal aquifer and the Tegelen fault zone (Table 4) and was later extended to also include the CLG doublet for the 2nd SHA. The model was updated after drilling and testing the CLG wells when a more detailed aquifer structure could be resolved (Table 1).
Table 4. Parameter for the hydraulic model used in the first two SHAs. The segmentation of the aquifer into two conductive and one intermediate, tight part reflects the identification of two major fluid loss zones at the CWG injection well. Both models differ in terms of the geometry, as the second model also includes the CLG-doublet and was extended in the Northern direction.
Simulation of thermal contraction stresses
Thermal contraction stresses were simulated on all known faults in the reservoir utilising parameters as listed in Table 2. An approach was chosen, where the cooled reservoir rock volume around the injection wells is approximated by (rectangular) elementary mechanical dislocation sources (Okada, Reference Okada1992). These sources represent the contraction of the cooled rock volume and facilitate the computation of resulting stresses at arbitrary points in the subsurface by superposition of contributions from the individual elementary sources. The approach significantly reduces the computational effort compared to a full-scale, thermo-mechanically coupled simulation procedure. At the same time, a high degree of accuracy can be achieved with this method. To demonstrate this, a generic doublet system in a horizontal, homogeneous aquifer is considered (model parameter according to Table 5). Thermal contraction stresses on a fault near the injection well were simulated with a thermo-mechanically coupled Finite Element model. The cooled reservoir rock (see Fig. 11) was also approximated by elementary sources, using two different spatial solutions. The results are compared in Fig. 12. While both approximations reproduce the stress pattern and magnitude of the Finite Element solution to first order, the high-resolution approximation yields a very good match also on a more detailed scale.
Table 5. Parameter used in the generic model of a geothermal doublet.
Fig. 11. Analytic solutions for the cooling front in a generic, homogeneous aquifer (after Schulz, Reference Schulz1987), displayed in a horizontal section at the centre of the reservoir layer. The cooled rock volume was simulated for different times after the start of the circulation according to the parameters in Table 5. The size of the elementary Okada sources is outlined for each time step by a red rectangle. Blue arrows indicate displacement directions outside the cooled zone, normalised for each time step. Black line denotes the trajectory of a fault located outside the cooled zone. Temperature decay with respect to the initial reservoir temperature according to the colour scale.
Fig. 12. Comparison of the Finite Element solution for thermal contraction stresses (Coulomb stress changes assuming a left lateral strike-slip failure mechanism) on a nearby fault with approximate solutions utilising Okada elementary sources. Approximation sources extend in the vertical direction over the complete reservoir layer. A segment of the cooled reservoir rock is approximated by three orthogonal, planar elementary sources, representing contraction of the volume in the three spatial directions. Approximation one utilises three elementary sources for the complete cooled rock volume with a lateral extension corresponding to the red rectangles as shown in Fig. 11. Approximation two consists of a larger number of elementary sources with a side length of 20 m, approximating the spatial shape of the cooled rock area shown in Fig. 11. Stress magnitudes according to the colour scale saturated at the minimum and maximum stresses for the Finite Element solution at each time step. The x-axis indicates distance along the strike direction of the fault.
Open access
