Validating induced seismicity forecast models - Induced Seismicity Test Bench
Eszter Kiraly-Proag, J. Douglas Zechar, Valentin Gischig, Stefan Wiemer, Dimitrios Karvounis, Joseph Doetsch
KKIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 1
Validating induced seismicity forecast models —Induced Seismicity Test Bench
Eszter Kir´aly-Proag , J. Douglas Zechar , Valentin Gischig , Stefan Wiemer ,Dimitrios Karvounis , and Joseph Doetsch An edited version of this paper was published by AGU. Copyright (2016) American Geophysi-cal Union. Kir´aly-Proag, E., J. D. Zechar, V. Gischig, S. Wiemer, D. Karvounis, and J. Doetsch(2016), Validating induced seismicity forecast models — Induced Seismicity Test Bench, J. Geo-phys. Res. Solid Earth, 121, doi:10.1002/2016JB013236. To view the published open abstract,go to http://dx.doi.org and enter the DOI.Corresponding author: E. Kir´aly-Proag, Swiss Seismological Service, ETH Zurich, NO H51.3,Sonneggstrasse 5, 8092 Zurich, Switzerland. ([email protected]) Swiss Seismological Service, ETH Zurich,Zurich, Switzerland. Swiss Competence Center for EnergyResearch (SCCER-SoE), ETH Zurich,Zurich, Switzerland
D R A F T D R A F T a r X i v : . [ s t a t . A P ] S e p - 2 KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
Key Points. ◦ We introduce a CSEP-based objective test bench for induced seismicityforecast models. ◦ We introduce a 3D smoothed seismicity model for induced earthquakes. ◦ We compare forecast models with different physical and statistical ele-ments on two EGS reservoirs.
Abstract.
Induced earthquakes often accompany fluid injection, and theseismic hazard they pose threatens various underground engineering projects.Models to monitor and control induced seismic hazard with traffic light sys-tems should be probabilistic, forward-looking, and updated as new data ar-rive. In this study, we propose an Induced Seismicity Test Bench to test andrank such models; this test bench can be used for model development, modelselection, and ensemble model building. We apply the test bench to data fromthe Basel 2006 and Soultz-sous-Forˆets 2004 geothermal stimulation projects,and we assess forecasts from two models: Shapiro and Smoothed Seismic-ity (SaSS) and Hydraulics and Seismics (HySei). These models incorporatea different mix of physics-based elements and stochastic representation of theinduced sequences. Our results show that neither model is fully superior tothe other. Generally, HySei forecasts the seismicity rate better after shut-in, but is only mediocre at forecasting the spatial distribution. On the otherhand, SaSS forecasts the spatial distribution better and gives better seismic-ity rate estimates before shut-in. The shut-in phase is a difficult moment forboth models in both reservoirs: the models tend to underpredict the seismic-ity rate around, and shortly after, shut-in.
D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 3
1. Introduction1.1. Induced seismic hazard
Seismicity caused by human activity, what is currently being called induced seismic-ity, is not a new phenomenon. Over the last several decades, workers have noted thatearthquakes are triggered by human activities including nuclear explosions [
Boucher et al. ,1969], fluid extraction [
Segall , 1989], fluid injection [
Seeber et al. , 2004;
Ellsworth , 2013],controlled filling of artificial reservoirs (e.g., Koyna, India) [
Gupta , 2002], and mining andexcavation [
McGarr , 1976]. But interest in induced seismicity has recently spiked, as hasthe rate of induced earthquakes in the central and eastern US [
Ellsworth , 2013;
Wein-garten et al. , 2015]. Here, it appears that fluid injections, primarily involving wastewater,are causing extensive seismic activity including events such as the 2011 m w . Kim , 2013], the 2011 m w . Horton ,2012], the 2011 m w . Keranen et al. , 2013], and the 2012 m w . Frohlich et al. , 2014].For modern deep geothermal energy projects, induced seismicity is a concern becausefluids must be injected to stimulate and enhance reservoir permeability, allowing the heatto be extracted. There are two recent examples in Switzerland: the Basel EGS experimentin 2006 [
H¨aring et al. , 2008] and the St. Gallen hydrothermal injection in 2013 [
Kraft et al. ,2013;
Edwards et al. , 2015;
Obermann et al. , 2015]. Both projects were canceled: Baselbecause of widely-felt seismic activity, and St. Gallen due to gas inflow, the low naturalfluid flow rate, and the high level of seismic activity during a short-term stimulation.These experiments demonstrated that project managers and operators have to be able tomanage induced seismic hazard and must strike a balance between reservoir creation (i.e.,
D R A F T D R A F T - 4
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS permeability enhancement, which is required for a geothermal system to be profitable)and induced seismicity. Induced seismicity during geothermal projects is a blessing anda curse: the spatial extent of micro-seismicity is a proxy for the size of the stimulatedreservoir, but felt and potentially-damaging earthquakes pose seismic risk to people andinfrastructure. Induced earthquakes in deep geothermal reservoirs are usually smallerthan m
3, but larger events ( > m
4) can occur, the largest so far being an m . Majer et al. , 2007]. Certainly, induced earthquakesfelt by the public may deter future geothermal projects. Despite the cancellations at Baseland St. Gallen, several geothermal projects in Switzerland are in development. As partof the Swiss national energy strategy, deep geothermal heat should supply 5 −
10% ofthe national baseload electricity [
Giardini , 2014]. One of the main obstacles to achievingthis goal is induced seismic hazard. To minimize induced seismic hazard, it is crucial notonly to monitor and analyze induced events, but also to develop a near-real-time toolfor making operational decisions. Such a hazard management scheme should be used toplan and operate reservoir stimulation so that large induced earthquakes are avoided [e.g.,
Bachmann et al. , 2011;
Mena et al. , 2013;
Goertz-Allmann and Wiemer , 2013].
Bommer et al. [2006] introduced a traffic light system to monitor and react to seismicactivity during geothermal reservoir stimulation. Like most traffic lights, this systemdistinguished three hazard levels, which were based on the size of events, observed peakground velocity, and public response. But the thresholds used to change the light werechosen subjectively, primarily by expert judgment [
Hirschberg et al. , 2015], and in practicethe system has resulted in operators taking action too late to avoid large events or a high
D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 5 seismicity rate. For example, in Basel the early induced earthquakes suggested that feltevents were likely, but the traffic light system failed to anticipate them [
H¨aring et al. , 2008].An improved hazard management scheme should be a dynamic, forward-looking systemthat incorporates real-time data and makes probabilistic forecasts of induced seismicityand its consequences. Such an Adaptive Traffic Light (ATL) system is composed of severalmodules (Figure 1):1. Collecting prior information, e.g., geological setting for hazard assessment and build-ing classifications for risk assessment (yellow in Figure 1). These data are essential to plana geothermal project and can address questions such as where to drill wells, the orien-tation of the local stress field, how to design reservoir creation plans, and the maximumpossible magnitude [
Gischig , 2015].2. Real-time data flow of hydraulic and seismic information (red in Figure 1). Theseare hydraulic data (e.g., injection flow rate and pressure measurements in the well) andseismic data that allow one to monitor reservoir creation, circulation, or other activitiesin the reservoir.3. Modeling and forecasting seismicity (orange in Figure 1). The key element in anATL system is seismicity forecasting. To forecast, we consider two periods: a learningperiod and a forecast period. During the learning period, seismic events are observedand analyzed according to their distribution in time and space. Then a calibrated modelforecasts the number, magnitude distribution, and spatial distribution of events in theforecast period.4. Ground motion models (gray in Figure 1). These models estimate the shaking thatan earthquake will cause and are based on properties of the earthquake source (e.g.,
D R A F T D R A F T - 6
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS its magnitude, style of faulting, and depth), wave propagation (distance to the earth-quake), and site response (type of rock, soil that can attenuate or amplify ground shaking).Ground Motion Prediction Equations [
Douglas et al. , 2013] and the Virtual EarthquakeApproach [
Denolle et al. , 2013, 2014] are examples of possible choices to estimate groundmotions.5. Combining models to account for epistemic uncertainties (green in Figure 1). Nosingle model captures all of the important features of seismicity. Model combination usingappropriate weights is one way to try to leverage each model’s best features.6. Calculating hazard and risk (brown in Figure 1). One can estimate the seismic haz-ard — the probability that some level of shaking will be exceeded — by combining groundmotion models and either synthetic catalogs generated by forecast models or individualscenario earthquakes. One can use this hazard to estimate the seismic risk: the potentialeconomic, social, and environmental consequences of seismicity.7. Guiding on-site decision-making processes (white in Figure 1). Based on hazard andrisk calculations, operators can make decisions concerning future stimulation strategiesand adjust flow rate accordingly.In this paper, we focus on the forecast models and the performance assessment modulesof the ATL system (delineated by a dashed gray line in Figure 1).
Induced seismicity models can be grouped into three classes [e.g.,
Gischig and Wiemer ,2013;
Gaucher et al. , 2015]: statistical, physics-based, and hybrid. In general, statisticalmodels for induced seismicity [e.g.,
Reasenberg and Jones , 1989;
Hainzl and Ogata , 2005;
Bachmann et al. , 2011;
Mena et al. , 2013] are conceptually and computationally simple
D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 7 and include aleatory uncertainty. But they do not explicitly account for the physicalprocesses governing induced seismicity (e.g., fluid flow in fractures, permeability changes,and stress interaction) and, until this study, they have not been used to forecast the spa-tial distribution of earthquakes. It is sometimes thought that statistical models, becausethey are primarily based on clustering, are limited in their ability to predict large eventsor make accurate long-term forecasts. In contrast, physics-based models [e.g.,
Olivellaet al. , 1994;
Bruel , 2005;
Kohl and M´egel , 2007;
Baisch et al. , 2010;
Rinaldi et al. , 2015;
McClure and Horne , 2012;
Wang and Ghassemi , 2012;
Karvounis and Wiemer , 2015;
Mignan , 2015] do consider underlying physical processes, and are hoped to perform bet-ter when operational conditions change, such as for the shut-in period, and for long-termforecasts. But the high computational expense of most physics-based models precludestheir use in near-real-time applications for the moment. Hybrid models are a compromisebetween physical models and statistical models. The goal of hybrid model developmentis to include some physical complexity and replace more complex physical considerationswith statistical methods or stochastic processes.
Mena et al. [2013] compared forecast models using the Basel dataset and found thatShapiro’s model [
Shapiro et al. , 2010] provided a good fit to the rate of induced earth-quakes. This model uses the seismogenic index, Σ, a parameter that describes the expectedseismic response of a given site. The seismogenic index is a function of the total injectedfluid volume and can be estimated from a short injection period or from the entire stim-ulation period; it also takes into account the b -value of the observed seismicity and thetotal injected volume. Using Σ, one can forecast the number of earthquakes in a givenmagnitude range and given period. Like most statistical models for induced seismicity D R A F T D R A F T - 8
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS [e.g.,
Bachmann et al. , 2011], Shapiro’s model does not make any predictive statementsabout the size or shape of the seismicity cloud. But it is crucial to monitor and anticipatethe shape and size of the seismic cloud during reservoir stimulation for two reasons. First,the extent of the seismicity cloud is used to estimate the volume of the stimulated reser-voir, which is crucial for energy production. Second, the spatial distribution of seismicityaffects hazard and risk analysis: many geothermal sites are located near settlements, mak-ing energy transportation cheap but posing a risk to infrastructure and people [
Edwardset al. , 2015]. Seismic risk strongly depends on geological settings (e.g., rock type underthe settlement), building vulnerability, and the depth of induced events. For instance, if a m w km below strong, new homes built on a rock site, almost all buildingswould remain intact, with only some slight damage. If an event of the same size occurs3 km below vulnerable houses built on a sedimentary basin, it is more likely that thehouses would be slightly damaged, and some houses may be moderately or even heavilydamaged [ Gr¨unthal , 1998]. Because the spatial distribution of induced seismicity is soimportant, any ATL system should be driven by 3D spatial forecasts.In this study, first we extend Shapiro’s model to produce 3D forecasts (SaSS model, i.e.,Shapiro and Smoothed Seismicity model). Then, we perform systematic statistical testson this model and on a hybrid model, in which seismicity is triggered by a numericallymodeled pressure diffusion (HySei model, i.e., Hydraulics and Seismicity model). To datethese are the only models in our institute, that are calibrated against real data, andsystematic re-calibration and testing can be carried out; moreover, they have a goodvariety of model features, which forecasts are worth evaluating and comparing. To dothis, we develop an Induced Seismicity Test Bench.
D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 9
Gerstenberger andRhoades , 2010;
Schorlemmer et al. , 2010;
Zechar et al. , 2010a;
Nanjo et al. , 2011;
Eberhardet al. , 2012;
Mignan et al. , 2013;
Taroni et al. , 2013;
Zechar et al. , 2013]. This supportcomes in the form of testing centers that CSEP operates; these centers allow modelersto check the consistency of their model with observations and to compare models. Wedescribe these activities in more detail in Subsection 3.2.The proposed Induced Seismicity Test Bench requires models to be tested, good qualityinduced seismicity datasets, and a robust statistical testing framework allowing objectivemodel evaluation. To test model consistency with observations and to rank models, we relyon pseudo-prospective forecasting, i.e., data that come from past stimulation experiments.Modelers calibrate their models using data recorded during a learning period and makeforecasts for a subsequent forecast period. Since observed data of the forecast periodsare already available, we can compare observed and forecast data after each recalibrationand test the consistency of the forecast in terms of seismicity rate, spatial distribution,and magnitude distribution. We can use statistical metrics such as the information gainper earthquake to compare model pairs and rank models according to their forecast skill[
Rhoades et al. , 2011]. Modelers should use the results of testing for further development,
D R A F T D R A F T - 10
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS creating a feedback between testing and modeling. The long-term goal is to develop anoperational ATL system to plan and conduct reservoir creation without a high rate ofseismicity or large events. A detailed flowchart of the Induced Seismicity Test Bench canbe found in the supplement (Figure S1).The Induced Seismicity Test Bench is a diagnostic tool: it can highlight which model ele-ments, be they physical or statistical, are essential for good forecasts, and why. This canin turn improve the models and our understanding of the underlying physical phenomena.In addition to using the test bench as a diagnostic tool, it can also be utilized on the flyto judge the performance of several models since the last forecast. The results can thenbe used for further improvement of the individual models and/or they can be applied toweight the models for the next forecast.In the next section, we briefly describe the data from two Enhanced Geothermal Sys-tems: the Basel 2006 experiment and the Soultz-sous-Forˆets 2004 stimulation. In Section3 we present two models, SaSS (Shapiro and Smoothed Seismicity) model and HySei(Hydraulics and Seismicity), which are calibrated on the datasets; and we also detail thetesting approach. We describe the testing results in section 4, discuss our findings insection 5, and conclude in section 6.
2. Data
The data we consider in this study come from the Soultz-sous-Forˆets 2004 and Basel2006 geothermal stimulations.The Basel geothermal site is located in northwestern Switzerland, at the southeastern partof the Upper Rhine Graben (Figure 2.a). The graben structure is an inactive extensionalrift system oriented N-S [
Zoback , 1992]. Here, the crystalline basement is covered by
D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 11 . km of sedimentary rock [ H¨aring et al. , 2008]. The well BASEL1 was drilled to a depthof 5 km between May and October 2006. In December 2006, after several hydraulic tests,the reservoir was hydraulically stimulated to enhance its permeability. The plan was tostimulate for 21 days, but after 6 days the injection was stopped due to intensive seismicity.In the year that followed, 3 additional events of m L > . H¨aring et al. , 2008].Based on the results of a subsequent risk study [
Baisch et al. , 2009;
Secanell et al. , 2009],the project was abandoned. After several years, the reservoir still has earthquakes, butthe seismicity rate is very low (1-3 earthquakes recorded per year) [
Deichmann et al. ,2014]. In this study, we use about 15 days of hydraulic [
H¨aring et al. , 2008] and seismicdata [
Dyer et al. , 2010] from the beginning of the stimulation (2006-12-02, 18:00), and wealso use the pre-stimulation injection test data.The Soultz-sous-Forˆets geothermal site is also located in the Upper Rhine Graben, betweenKutzenhausen and Soultz-sous-Forˆets, about 70 km north of Strasbourg (Alsace, France;inset in Figure 2). The geothermal gradient is about 100 ◦ C/km within the 1 . km thicksedimentary cover over a granitic basement [ Evans et al. , 2012]. This abnormally highgeothermal gradient is related to deep hydrothermal convection cells in the fracturedbasement [
G´erard et al. , 2006]. The geothermal project here started in the early 1980s andfour wells have been drilled into two reservoirs: one at about 3 . km depth (GPK1, GPK2wells) and another at about 4 . km (GPK2, GPK3, GPK4 wells). Several stimulations andcirculation tests were carried out [ G´erard et al. , 2006;
Cal`o et al. , 2014;
Genter et al. , 2012].Energy production started in 2008 [
Genter et al. , 2010]. In this study, we use hydraulic andseismic data of the pre-stimulation and stimulation of September 2004 (Figure 2.b,
Dyer [2005]). Local magnitudes were corrected by using the scaling relationship by
Douglas
D R A F T D R A F T - 12
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS et al. [2013]. Note that the seismograms in this data set are clipped, causing saturationof the magnitudes at 1 .
8; that is, no event has m w > .
3. Models and testing3.1. The Shapiro and Smoothed Seismicity (SaSS) model
The SaSS model is computationally simple and based on the seismogenic index, Σ[
Shapiro et al. , 2010]; we distribute the earthquakes expected by Σ in 3D by smoothingseismicity in space. Shapiro’s model, which describes the rate of induced seismicity duringstimulation, is defined as: log ( N m ( t )) = log ( Q c ( t )) − bm − Σ (1)where N m ( t ) indicates the number of induced events above magnitude m up until time t , Q c ( t ) denotes the cumulative injected volume of fluid at time t , b is Gutenberg-Richter b -value of the observed seismicity, and m is the magnitude above which all events areexpected to be reliably recorded (often called the magnitude of completeness).To forecast the number of events in the forecast period, we estimate Σ and b from thelearning period, and we predict the total volume that will be injected by the end of theforecast period. Kir´aly et al. [2014] compared four deep geothermal datasets and foundthat in some cases b and Σ are not constant during and after stimulation; thus, we re-estimate them at the end of each learning period, every six hours. To predict Q c ( t ) at theend of a forecast period, we assume that the injection flow during the forecast period willfollow the previously-planned strategy. Eq. 1 describes the rate of induced seismicity onlyduring stimulation [ Shapiro et al. , 2010]. As soon as the stimulation stops (the momentof well shut-in), the rate of induced earthquakes is expected to decay; the SaSS model
D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 13 assumes the decay follows the equation of
Langenbruch and Shapiro [2010] (using theoriginal notation for consistency): R b (cid:18) tt (cid:19) = R a (cid:18) tt (cid:19) p (2)where R b is the post-stimulation seismicity rate at time t (since the beginning of thestimulation), t is the length of the stimulation period before shut-in, R a denotes the av-erage seismicity rate during stimulation, and p controls how quickly the rate decays. Forsubsequent forecast time windows (i.e., 6-hour time bins of the forecast period, FTWs),the majority of parameters are calibrated on the corresponding learning period, but Q c and R b are recalculated for each time window. If the learning period ends in the stimula-tion period but some FTWs expand to the post-stimulation, the estimation of parameter p is not possible, thus we use a generic value: p = 2. Also, if p is estimated to be smallerthan 2 we set the value to 2, following the value that is proposed by Langenbruch andShapiro [2010] for an early post-injection period. Detailed flowchart of number componentcan be found in the supplement (Figure S2). As in CSEP experiments and suggested by
Shapiro et al. [2010], the number of events in each forecast period is assumed to follow aPoisson distribution and the numbers obtained by using Eq. 1 and 2 are Poisson expectedvalues; error bars in all subsequent figures indicate the 95% Poisson confidence interval.To model the 3D spatial distribution of induced earthquakes, we added a spatial compo-nent to the model by smoothing the seismicity observed during the learning period (Fig-ure 3.A). Several studies, including the Regional Earthquake Likelihood Models (RELM)experiment [
Schorlemmer et al. , 2010;
Zechar et al. , 2013] have shown that smoothedseismicity models are effective at forecasting the spatial distribution of tectonic earth-
D R A F T D R A F T - 14
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS quakes. To construct a smoothed seismicity model in two dimensions, one applies a two-dimensional smoothing kernel to each past event [e.g.,
Helmstetter et al. , 2007], calculatesthe contribution of smoothed earthquakes on a given grid, then sums contributions of allobserved earthquakes. To create a probability density function (PDF, i.e., earthquakespatial probability map), one normalizes the smoothed seismicity map so its sum is unity.We extend the 2D Gaussian smoothed seismicity model of
Zechar and Jordan [2010]to 3D. For each forecast period, we smooth all prior events, where the contribution of anearthquake to a given voxel (i.e., volume element) is K ( x e , y e , z e , x , x , y , y , z , z ) = 18 (cid:34) erf (cid:18) x − x e σ √ (cid:19) − erf (cid:18) x − x e σ √ (cid:19)(cid:35) × (cid:34) erf (cid:18) y − y e σ √ (cid:19) − erf (cid:18) y − y e σ √ (cid:19)(cid:35) × (cid:34) erf (cid:18) z − z e σ √ (cid:19) − erf (cid:18) z − z e σ √ (cid:19)(cid:35) (3)where x e , y e and z e denote the location of the given earthquake, x , x , y , y , z and z are the points that define the edges of the voxel, and σ , σ and σ are bandwidths ofthe 3D Gaussian kernels in EW, NS and vertical directions, respectively. To make a goodsmoothed seismicity forecast, we need good bandwidths; we optimize these by dividingdata from the current learning period into a training set and a validation set (Figure 3.C).The length of the training and validation sets depend on the length of the forecast periodand the learning period. If the length of the forecast period is more than half the lengthof the learning period, the training and validation sets are each one-half of the learningperiod. Otherwise, the length of the validation set is equal to the length of the forecastperiod. We search for the bandwidth combination that, when used to smooth the training D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 15 set, best forecasts the seismicity of the validation set. To avoid ’surprises,’ i.e., eventsoccurring where the model would not expect any events, we distribute a certain fractionof the PDF over all voxels (i.e., surprise factor), following the idea of
Kagan and Jackson [2000]. We analyze the performance of 1000 combinations of bandwidths and surprisefactors using the training and validation set of the learning period. The PDF is updatedfor each new learning/forecast period. Since the PDF is based on the learning period, thismodel assumes that earthquake locations in the forecast period will not be very differentfrom the seismicity observed so far.Smoothed induced seismicity models must differ from their tectonic counterparts in atleast one aspect: induced models should capture the propagation of the seismicity frontafter shut-in. In particular, due to pore pressure diffusion, induced seismic activity tendsto decrease in the vicinity of the injection well and to concentrate at the boundaries ofthe reservoir. We attempt to model this time-dependent effect by applying exponentialtemporal weighting: the most recent event receives a maximum weight (one), and earlierevents get smaller weights according to their origin time. This is analogous to the ex-ponential smoothing approach commonly used in time series forecasting [
Goodwin , 2010]and is also connected to the Omori-Utsu relation describing aftershock decay rate [
Zhuanget al. , 2012].The forecast magnitude distribution is the Gutenberg-Richter distribution [
Gutenbergand Richter , 1944] with the b -value estimated from the learning period. The HySei model developed by
Gischig and Wiemer [2013] describes seismicity triggeredby pressure diffusion with irreversible permeability enhancement. The biggest advantage
D R A F T D R A F T - 16
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS of the model is that it quantifies permeability enhancement by calibrating flow rate andwellhead pressure against observations. The HySei model consists of two main parts:hydraulic inversion and seismicity modeling. The aim of inverting hydraulic observationsis to reconstruct the pressure evolution in the reservoir. We seek the best hydraulicparameters to match the observed well-head pressure with a one-dimensional radial flowmodel. We use a finite difference method in a circle of 1200 m radius distributed on 3000nodes, and 1-minute resolution in time. During the pre-stimulation test injection, we solvethe diffusion equation (Eq. 4) with constant permeability ( κ = κ ). During stimulationthe governing equations are the diffusion equation (Eq. 4) with irreversible changingpermeability (Eq. 5) due to increasing pressure that exceeds some threshold (Eq. 6): ρS ∂p∂t = ∇ (cid:16) κρµ ∇ p (cid:17) + q m (4) κ = κ ( u + 1) (5) ∂u∂t = C u H pt (cid:16) ∂p∂t (cid:17) H u ( u t − u ) H p ( p − p t ) (6)where ρ is fluid density, S is the specific storage coefficient, κ is permeability that variesduring the stimulation, µ is fluid viscosity, and q m is a mass source; κ is the initialpermeability before the stimulation, u is stimulation factor (i.e., the overall permeabilityenhancement of the reservoir); C u is stimulation velocity, a constant that scales the rateat which permeability changes, u t is maximum stimulation factor, and p t is thresholdpressure, H pt is a Heaviside function, it is one if pressure increases, zero otherwise, H p and H u are Heaviside functions for pressure and stimulation factor. These are smoothed D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 17 to avoid a singularity and resulting numerical instability. Permeability starts to increase ifpressure reaches p t . If pressure further increases, the permeability of the reservoir increasesuntil it reaches u t . Note that a reversible component of permeability change representingthe compliant response fracture to pressurization [e.g., Rutqvist and Stephansson , 2003]has not been included in this version of the model.In the seismicity model, randomly-placed potential nucleation points are triggered bythe radial symmetric pressure evolution following the Mohr-Coulomb failure criterion.They have no spatial extent, but differential stress ( σ − σ ) is defined at the seed point.Local b -values are determined at the seed points following a linear relationship betweendifferential stress and b -value: b max and b min parameters are b -values at minimum and max-imum values of differential stress, respectively. When a seed point is triggered, a randommagnitude is drawn from the magnitude distribution with the local b -value. Additionalfree parameters are the scaling factor F s (the ratio between the number of synthetic andobserved events), the stress drop coefficient dτ (the change of stress conditions after aseed has been triggered), and a criticality threshold dµ , which accounts for the fact thatseed points cannot be too close to the failure limit.For this study, we parallelized parts of the code and extended the model to 3D (Figure3.B) by adding an off-fault component to the originally 2D seismicity model. Assumingthat the seismicity is generated on the current main fault, we determine the principalcomponents of the current seismicity cloud and use the empirical distribution of the seis-micity along the smallest axis to define off-fault coordinates of the synthetic events. Adetailed flowchart of the HySei model can be found in the supplement (Figure S3).To represent the spatial differences of the two models, Figure 4 shows cross sections of D R A F T D R A F T - 18
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS the 3D PDFs of SaSS (upper line) and HySei (bottom line) at the moment and locationof the biggest event ( m w = 3 . To assess a single model, we check if its forecasts are consistent with the observations[
Zechar et al. , 2010b], asking the question: might the observations have been generatedby this model? One way we do this is to check if the number of observed earthquakesfalls within the 95% confidence interval of the forecast. If so, the model passed theNumber-test. In a similar way, we examine if the magnitude distribution of all forecasts isconsistent with the observations (Magnitude-test). To test the spatial component (Space-test) [
Zechar et al. , 2010b;
Rhoades et al. , 2011], we use a testing grid of 4 km × km × km centered on the well tip and divided into 200 m × m × m voxels. After normalizingthe forecasts so that the number of forecast events matches the number of observed events,we calculate the log-likelihood (LL) of the observation in each voxel. Summing these valuesgives a joint LL for a specific experiment. The higher the joint LL values are the betterthe forecast [ Zechar et al. , 2010b;
Rhoades et al. , 2011].To check if the forecast is consistent with the observed seismicity of the forecast period,we simulate 1000 catalogs from the forecast, and find the 5 th percentile of the LL valuesfor the simulated catalogs. If the LL for the current observation is higher than the 5 th percentile the forecast passed the Space-test — the observed seismicity could have beengenerated by the model. Both models consider the earthquake distribution Poissonian,thus LL values are calculated as follows: D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 19 L ( A ) = n (cid:88) i =1 (cid:104) k i × log (cid:0) λ A i (cid:1) − λ A i − log (cid:0) k i ! (cid:1)(cid:105) (7)where L ( A ) is the Poisson joint LL of forecast A, n is the number of voxels, k i is thenumber of earthquakes observed in the ith voxel, and λ A i is the forecast seismicity ratein the i th voxel of forecast A.To compare two models, one can directly compare individual LL values of the modelseither for model components (i.e. event numbers, magnitudes or the spatial component)separately or for the entire model. These measures give information about the modelperformance not only against data but against other models. Here we would like toemphasize that LL values consider the whole model space. In other words, it reflectsthe performance of not only the temporal/magnitude/spatial bins that host at least oneearthquake but also the empty ones answering the question: what is the probability tohave zero earthquake in the given temporal/magnitude/spatial bin?One can also calculate the information gain of one model with respect to another formodel comparisons. This measure emphasizes the non-empty bins by comparing theforecast seismicity rates of model A with that of model B in the voxels where earthquakesoccurred. The following formula gives I i , the information gain of model A over model B for an earthquake occurring in the i th voxel [ Rhoades et al. , 2011]: I i = − N A + N B N + ln (cid:32) λ A i λ B i (cid:33) (8)where N is number of observed events, λ A i and λ B i denote forecast seismicity rate inthe i th voxel of model A and B , respectively, N A and N B are the total forecast numberof events in model A and B , respectively. The first term of the right hand side is a D R A F T D R A F T - 20
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS penalty concerning the number of events under each model. We seek to know if onemodel is better than the other, in other words, if the expected value of the informationgain population differs from zero. One can also estimate how much better or worsemodel A relative to model B (i.e., average information gain) by finding an appropriateestimator. Exponentiating the average information gain yields the average probabilitygain of model A with respect to model B . Additionally, 95% confidence interval of theestimated expected value can be calculated to determine if model A is significantly betteror worse than model B : if the confidence interval contains zero, the difference betweenthe models is not statistically significant at the 5% significance level.Several techniques are possible to compute the average information gain. Rhoades et al. [2011] suggested to take the arithmetic mean of the information gain distribution asthe expected value of the population, based on Students t-distribution [
Student , 1908].We refer to this method as ’Classical mean’. This estimator is best if the populationfollows a normal distribution. Plotting the distribution of information gains (that is, forindividual earthquakes) for SaSS relative to HySei as a function of time and in a quantile-quantile plot (Figure S4) suggests that the information gains are not normally distributed.One possible way to solve this problem is to seek an estimator that can tackle outlierssystematically. This can be done by manual data screening and removal of outliers,but it can be impractical due to the large number of data points and possible masking(i.e., large outliers can hide smaller ones). To overcome these problems, we use robuststatistics to automatically detect and downweight outliers [
Ruckstuhl , 2014]. We refer tothis method as ’Robust mean’. To calculate the expected value of the information gaindistribution, we compute a weighted mean where the influence of the outliers is reduced.
D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 21
In particular, we use the Huber M-estimator, implemented as mlochuber in the LIBRAmatlab package [
Verboven and Hubert , 2005]. By using the Huber M-estimator, we avoidthe problem that a few earthquakes dominate the estimate of the average informationgain. We also explore a non-parametric method: generate 1000 bootstrap samples of theobserved information gains (i.e., we sample with replacement) and find the arithmeticaverage and 2 .
5% and 97 .
5% percentiles, thus obtaining a ”Bootstrap mean” and thecorresponding 95% confidence interval. Using the same bootstrap samples we also find’Bootstrap median’. We show a comparison of these methods in the next section.
4. Results4.1. Consistency tests
Figure 5 shows four snapshots of forecast and observed seismicity rates for both datasets.The top row shows the corresponding hydraulic data (injection rate and well-head pres-sure) to provide time reference for the forecasts. Blue, red, green, and purple verticallines indicate the end of the different learning periods: corresponding shaded areas showforecasts of SaSS model (middle row) and HySei model (bottom row) with 95% Poisso-nian confidence intervals. In case of Basel 2006, both models seriously overpredicts theseismicity rate for LP1 (blue learning period that ends at day 1 . . .
25) andLP4 (purple learning period that ends at day 9 . . D R A F T D R A F T - 22
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS underpredicts (after LP2 the learning period that ends at day 3 . .
5, respectively). HySei performs well in most of the cases (after LP2, LP3 andLP4), except after LP1. In this case, the model expects higher pressure in response tothe injection peaks between day 2 −
3, which results in overprediction of the sesimicityrate. This might be due to the fact that a reversible component of permeability change,possibly arising from fracture compliance, is not included in this version of the model.To show forecasts corresponding to all learning periods, we use a matrix representationwhere colors indicate the goodness of the forecast (Figure 6): yellow means a perfect fore-cast; red and blue mean under- or overprediction, respectively. Downward- and upward-pointing triangles denote moments when the observed seismicity rate falls out of the 95%confidence intervals due to serious under- or overprediction, respectively. To avoid overlapof the forecast periods, we represent the 3-day forecast period vertically: the end of thelearning period is indicated on the horizontal axis, time during the 3-day forecast period isindicated on the vertical axis with subsequent 6-hour FTWs. Time in the forecast periodincreases from bottom to top. The top row of Figure 6 shows the observed seismicity ratefor both datasets, middle and bottom rows show a comparison of observed seismicity rateswith forecasts from SaSS and HySei, respectively. In Basel, both models mainly overes-timate the number of observed earthquakes during the initial stimulation period. Whenthe injection rate was decreased and at shut-in, both models have difficulties forecastingthe right number of earthquakes: they severely underpredict the observed seismicity rate.The SaSS model overpredicts for the post-stimulation period, whereas HySei seems to findgood estimates most of the time for later periods (with the exception of three time win-
D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 23 dows). In Soultz-sous-Forˆets 2004, the SaSS model mainly forecasts well or overestimatesthe number of earthquakes during stimulation. The forecast period corresponding to thelearning period of day 3 . . . . m > .
8. Both models forecast the magnitude distribution of micro-seismicevents well, meaning that observed seismicity follows the Gutenberg-Richter relation inalmost all cases. Nevertheless, the probability of the biggest event of the Basel 2006project is very small in both models (insets in Figure 7b-c). The truncated magnitudes inSoultz-sous-Forˆets 2004 preclude us from considering the probability of the largest eventin this data set, because we have no good estimate for the magnitude of the largest event.
D R A F T D R A F T - 24
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
We investigate the spatial component of the models by dividing the joint LL by thenumber of observed events (LL/Eqk) in Figure 8. We decided to normalize due to thefact that LL values are correlated with the number of earthquakes in a FTW. We use thesame matrix representation as we introduced for the number component: end of learningperiods are indicated on the horizontal axis, FTWs on the vertical axis. Yellow indicatesbetter results than red, the higher the LL value, the better the forecast is. Crossesrepresent moments when the model does not pass the Space-test. Gray squares denotemoments when no earthquake occurred. Gray dotted line marks the shut-in moment. It isclear that SaSS passes the Space-test more often than HySei does, especially after shut-infor both datasets. Additionally, SaSS’s LL values are higher than that of HySei indicatingthat smoothed seismicity outperforms the simple geometry of HySei’s forecasts.
To be able to compare the two models we calculate LL from the absolute values ofthe Number- and Magnitude-test by answering the same question we addressed in caseof the spatial component: what is the probability of the observation given the modelforecast? We calculate LL values for all FTWs of all model components (Figure S5-S6).Figure 9 gives an overview of differences between the model LLs. Green shows when SaSSperforms better than HySei, pink shows when HySei is better than SaSS, white indicatesthat the models forecast similarly. The magnitude component is exceptional in this figure,because we do not test the consistency of the forecast and observations in incrementalFTWs, rather the cumulative distribution. For instance, in case a 3-day magnitude testwe take all events occurred in the forecast period from the end of the learning period until
D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 25 the end of day 3. This allows a more stable distribution of the observed events that canbe tested against a power law.These results clearly confirm that the magnitude component is very similar in the mod-els, which is not surprising since both models use the Gutenberg-Richter relation. Thedifferences lay in the number and spatial components. In terms of number, SaSS performsbetter in several moments during the stimulation and in the early post-stimulation periodin Basel. HySei gives better results close to the shut-in and generally after the stimu-lation, especially at later moments of the experiment. The green color in most FTWsof the spatial component reveals that SaSS holds the better spatial component, which isemphasised towards the end of the experiment.To compare the entire model performance, we merge all components and calculate LLnormalized by the number of earthquakes occurred in the given FTW. Figure 10 detailsthe sum of LL/Eqk values of the individual FTWs for 6-, 24-, 48-, and 72-hour forecastperiods. Three regimes can be observed in the case of Basel 2006: • regime A : when models perform similarly well • regime B : when SaSS model is better than HySei • regime C : when HySei overcomes SaSS, especially for the longer forecast periods.Comparing these results to the performance of individual model components, it is clearthat the regimes are determined by the interplay of the number and spatial components.Both components of both models perform similarly in regime A , which results in similaroverall performance. Around the shut-in, even if HySei gives better number forecasts fora short period, SaSS can compensate with its spatial component and it overcomes HySeialso with its number component by the end of regime B , which results in a better overall D R A F T D R A F T - 26
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS performance of SaSS for this period. As the number of events drastically decreases relativeto previous periods in regime C , it seems that HySei’s more precise number forecastscompensate against SaSS’s better spatial forecasts giving better overall LL values. In thecase of Soultz-sous-Forˆets 2004, only two of the three regimes are present: regime B fromthe beginning of the experiment about 1 . C for the rest of the experiment. In the first part ofregime B , the slightly better spatial component of SaSS compensates the generally betternumber component of HySei giving marginally better results to SaSS. From the shut-inuntil the end of regime B , the spatial component of SaSS is clearly better together withthe fact that HySei’s number component is less dominant than previously. This results ina drop of overall LL. The decrease of number of induced earthquakes (regime C ) highlightsagain that HySei’s number component overcomes SaSS’s better spatial component.Summarizing the model comparison based on LL: SaSS obtains better results in spacegenerally, in terms of seismicity rate in some moments of the stimulation, and also theentire SaSS model gives better results until a certain point after shut-in (regime B ) for bothdatasets; HySei outperforms SaSS in seismicity rate forecast in the post-stimulation periodand also the overall LL values of HySei in the late post-stimulation period, especially forlonger forecast periods.Figure 11 presents the results of all 6-hour information gains from the beginning until theend of the experiment for both datasets. Solid black lines indicate the empirical probabilitydensities of the information gains, dotted gray lines denote normal distributions, where theexpected values and standard deviations are estimated from the corresponding empiricaldistributions. To use the classical method to determine the average information gain, D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 27 the population should be normally distributed. This is not the case, which is why weinvestigate four methods to calculate the average information gain: classical mean, robustmean, bootstrap mean, and bootstrap median corresponding to red, green, orange, andbrown, respectively. Insets show the estimated average values with their uncertainties.For both datasets medians and robust means are closer to the the clear peaks of thepopulations, whereas classical and bootstrap mean values are shifted and have wider95% confidence intervals. In the case of the Basel 2006 data, interpretation of modelperformance depends on the choice of the estimator: for robust mean and bootstrapmedian HySei performs significantly better than SaSS, for classical and bootstrap meanexactly the opposite. This emphasizes that we should be cautious about information gaininterpretations. In our opinion, in case of information gain calculations, (1) it is necessaryto check the distribution of the observed information gains, (2) it is recommended to useseveral estimators to have a clearer view of the possible average information gain values,and (3) to interpret the results carefully. An overview of average information gain for6-, 24-, 48, and 72-hour forecast periods with all four estimators can be found in thesupplement (Figure S7-S10).
5. Discussion
Predictive models of induced earthquakes can help reduce seismic hazard and risk duringreservoir stimulations. Although many models are being developed, most are presentedin a context that is descriptive, not predictive: they are tuned using the entire data set,and so their ability to forecast is not checked. In this study, we propose a test bench toobjectively evaluate various induced seismicity models. We bring to the test bench twomodels used to forecast two datasets. We demonstrate that such a test bench can quantify
D R A F T D R A F T - 28
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS the forecast skill of different models. The results can give guidance how to merge models.One possible way to combine models is weighting models by their past performance. Thetest bench can provide detailed information about the performance of the tested modelsthat can be converted to probabilistic weights. Weighted average models has the potentialto merge the best forecasting features of the tested models and can give important inputfor real-time forecasting and hazard assessment. The test bench can also highlight modelfeatures to be improved, e.g., because the model performs badly at forecasting one of thekey parameters (i.e., event number, magnitude distribution, or spatial distribution) orduring certain moments (e.g., during stimulation, at shut-in, or after shut-in).Our test bench showed that both tested models are limited to accurately forecast therate of induced earthquakes. The forecasts are particularly bad around shut-in. Duringstimulation and shortly after shut-in, we observe first a slight overprediction and then asevere underprediction as the injection rate decreases and stops. In the post-injection pe-riod, SaSS overpredicts the number of events (except the moment when model parametersare not well calibrated due to the very short post-injection period).As suggested by
Langenbruch and Shapiro [2010], we use a generic value of 2 for pa-rameter p when parameter estimation is not possible, and the same generic value is usedif calculated ones are lower than 2. In Basel, we observed that calculated values of p arealways smaller than 2. This means that we always apply a decay with p = 2, which resultsin faster decay than the data of learning period would suggest. Nevertheless, all modeleddecays are slower than the observed seismicity decay, indicated by massive overpredictionsin the post-stimulation periods. In contrast, for Soultz-sous-Forˆets estimated values of p are always higher than 2 allowing good forecasts at the beginning of the post-stimulation D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 29 period but the decreasing tendency of the values of p results in overpredictions for laterforecast periods. These results suggest that forecasting the post-injection seismicity isdifficult and the current post-injection seismicity decay law is not appropriate in an op-erational forecasting environment.The spatial forecasts of the SaSS model gave generally good results. But these forecastsare limited by the fact that they are based on the current learning period. The model cangive good forecasts when the seismicity is nearly stationary, i.e., new earthquakes occurwhere previous ones occurred. But this is often not the case in induced seismicity relatedto geothermal reservoir creation, where seismicity propagates with the pressure front.In future work, to incorporate diffusion-like propagation of the seismicity, we imagine astep-by-step spatial forecast for each FTW of the forecast period. One could simulatethousands of synthetic catalogs for the first FTW based on the learning period. Forecastsof FTWs are based on the PDF calculated from the synthetic catalogs of the previousFTWs. Temporal weighting (exponential or some other temporal weighting) of generatedearthquakes can help to simulate the migration of the seismicity cloud.One might also improve induced seismicity forecasting by considering Coulomb stresschanges, which has been shown to a good descriptive model of tectonic seismicity [ Steacyet al. , 2005] and has been considered in the induced seismicity context:
Orlecka-Sikora [2010] suggested that static stress transfer can have an accelerating impact on mining-induced seismicity, and
Schoenball et al. [2012] concluded that static stress change doesnot play an important role during stimulation but might help to trigger after shut-in in theSoultz-sous-Forˆets reservoir. Moreover,
Catalli et al. [2013] found that 75% of the analyzedinduced earthquakes (based on
Deichmann and Ernst [2009]) in Basel occurred in regions
D R A F T D R A F T - 30
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS of increased Coulomb stress, where failure is thought to be encouraged. Unfortunately,prospective tests of the Coulomb stress hypothesis are difficult because one needs accurate,real-time estimates of hypocenter, magnitude, and focal mechanism, and one also needssome a priori knowledge on fault orientations in the reservoir.Additional model improvements may relate to the statistical description of earthquakedistributions. In the testing framework and also in all CSEP experiments, earthquakeoccurrence is considered as a Poissonian process [
Eberhard et al. , 2012]; LL and confidenceinterval computations are based on that assumption. The Poissonian assumption is notcompletely fulfilled, because earthquakes are not independent, neither in time nor inspace.
Eberhard et al. [2012] reported that Poissonian distribution was not supported bythe seismic data; others [e.g.,
Kagan , 2010;
Lombardi and Marzocchi , 2010] have previouslyshown the same observation in different regions and magnitude ranges. Failures of modelforecasts might stem from the Poissonian assumption beside the fact that the model doesnot incorporate the necessary physical processes. Modeling earthquake occurrence as aPoissonian process is thus not ideal and improvements are subject of further investigations.It is necessary to emphasize that all tests are highly dependent on the observed catalog.Thus, it is extremely important to detect events and to determine good origin times,magnitudes and precise locations. For the moment, it is still a challenge, especially innear real-time.Our analysis further revealed that forecasting the rate and magnitude distributionsaround shut-in also remains a difficult question: the models often underpredict duringthis period and do not represent the magnitude distribution well. Presumably, this prob-lem is not specific to the data we considered here because in several other projects the
D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 31 biggest event occurred after shut-in [
Baisch et al. , 2006;
Asanuma et al. , 2005]. Focusingon shut-in and the events that follow,
Barth et al. [2013] showed theoretically and alsoconfirmed with the analysis of the data from Soultz-sous-Forˆets 2000 that probabilityof exceeding a certain magnitude can be higher after shut-in than it would have beenfor on-going injection.
Segall and Lu [2015] proposed a descriptive model that includescomplete poroelastic coupling — changes in pore pressure induce stresses, and changesin mean normal stress induce changes in pore pressure — and concluded that an abruptshut-in can produce sharp increase in the seismicity rate. Post shut-in peaks of the seis-micity rate result from the rapid change in stress before the pore pressure can be relieved.Concerning post shut-in magnitudes,
Segall and Lu [2015] claimed that larger events areabsent at short injection times but as injection proceeds the probability of larger earth-quakes increases, thus larger events occurring post shut-in are not unexpected. Anotherexplanation for large post-stimulation events came from
McClure [2015]: simulation withthe three-dimensional version of CFRAC [
McClure , 2012] revealed that post-stimulationseismic events can be caused by backflow from dead-end fractures into fractures that hostthe largest event. He proposed that pumping of fluid to the surface immediately aftershut-in could mitigate this effect and reduce post-stimulation seismic activity. The in-ferences made from these descriptive models ought to be used in future work to improvepredictive models such as those considered in this study.
6. Conclusions
Forward-looking, near-real-time warning systems can help avoid large induced earth-quakes and keep micro-seismicity at a tolerable level during and after project operations.The Induced Seismicity Test Bench can be used to test the core of such a warning system,
D R A F T D R A F T - 32
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS an Adaptive Traffic Light system. Here, we tested, compared and ranked the performanceof the SaSS and the HySei models.To say which of these models performs best is not straightforward. In terms of magni-tude, both models forecast micro-seismicity fairly well, but none of them is able to forecastthe biggest m w . D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 33 test their models for consistency, and compare model performance: we believe this is themost efficient way to reduce induced seismic hazard.
Acknowledgments.
We acknowledge the GEOTHERM, GEOTHERM-2 and GEIS-ERS projects for financial support to develop fundamental ideas concerning the AdaptiveTraffic Light System and providing the stimulation data of Soultz-sous-Forˆets 2004. Theauthors would like to thank EEIG Heat Mining for permission to publish the data. Ac-knowledgement is also due to the numerous agencies which have supported the Soultzproject over the years including the European Union, ADEME of France, BMU of Ger-many, SER and SFOE of Switzerland, and the EEIG ’Exploitation Mini`ere de la Chaleur’consortium. Access to the data is provided by contacting the authors. We thank ArnaudMignan, Antonio Pio Rinaldi and Eduard Kissling for their valuable comments on anearlier version of the manuscript. We also thank Yehuda Ben-Zion as editor, the associateeditor, Carsten Dinske and three anonymous reviewers for their comments and sugges-tions. E.K.-P. acknowledges the GEOTHERM-2 project for financing her PhD. This workhas been partially completed within the Swiss Competence Center on Energy Research- Supply of Electricity, with the support of the Swiss Commission for Technology andInnovation.
References
Asanuma, H., H. Nozaki, H. Niitsuma, and D. Wyborn (2005), Interpretation of mi-croseismic events with larger magnitude collected at Cooper Basin, Australia,
GRCTransactions , , 87–91. D R A F T D R A F T - 34
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
Bachmann, C. E., S. Wiemer, J. Woessner, and S. Hainzl (2011), Statistical analysis ofthe induced Basel 2006 earthquake sequence: introducing a probability-based moni-toring approach for Enhanced Geothermal Systems,
Geophysical Journal International , (2), 793–807, doi:10.1111/j.1365-246X.2011.05068.x.Baisch, S., R. Weidler, R. V¨or¨os, D. Wyborn, and L. de Graaf (2006), Induced seis-micity during the stimulation of a geothermal HFR reservoir in the Cooper Basin,Australia, Bulletin of the Seismological Society of America , (6), 2242–2256, doi:10.1785/0120050255.Baisch, S., R. V¨or¨os, R. Weidler, and D. Wyborn (2009), Investigation of fault mech-anisms during geothermal reservoir stimulation experiments in the Cooper Basin,Australia, Bulletin of the Seismological Society of America , (1), 148–158, doi:10.1785/0120080055.Baisch, S., R. V¨or¨os, E. Rothert, H. Stang, R. Jung, and R. Schellschmidt (2010),A numerical model for fluid injection induced seismicity at Soultz-sous-Forˆets, In-ternational Journal of Rock Mechanics and Mining Sciences , (3), 405–413, doi:10.1016/j.ijrmms.2009.10.001.Barth, A., F. Wenzel, and C. Langenbruch (2013), Probability of earthquake occurrenceand magnitude estimation in the post shut-in phase of geothermal projects, Journal ofSeismology , (1), 5–11, doi:10.1007/s10950-011-9260-9.Bommer, J. J., S. Oates, J. M. Cepeda, C. Lindholm, J. Bird, R. Torres, G. Mar-roqu´ın, and J. Rivas (2006), Control of hazard due to seismicity induced by ahot fractured rock geothermal project, Engineering Geology , (4), 287–306, doi:10.1016/j.enggeo.2005.11.002. D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 35
Boucher, G., A. Ryall, and A. E. Jones (1969), Earthquakes associated with under-ground nuclear explosions,
Journal of Geophysical Research , (15), 3808–3820, doi:10.1029/JB074i015p03808.Bruel, D. (2005), Using the migration of induced micro-seismicity as a constraint for HDRreservoir modelling, in Thirtieth Workshop on Geothermal Reservoir Engineering , pp.1–7.Cal`o, M., C. Dorbath, and M. Frogneux (2014), Injection tests at the EGS reservoir ofSoultz-sous-Forˆets. Seismic response of the GPK4 stimulations,
Geothermics , , 50–58,doi:10.1016/j.geothermics.2013.10.007.Catalli, F., M. A. Meier, and S. Wiemer (2013), The role of Coulomb stress changesfor injection-induced seismicity: The Basel enhanced geothermal system, GeophysicalResearch Letters , (1), 72–77, doi:10.1029/2012GL054147.Deichmann, N., and J. Ernst (2009), Earthquake focal mechanisms of the induced seismic-ity in 2006 and 2007 below Basel (Switzerland), Swiss Journal of Geosciences , (3),457–466, doi:10.1007/s00015-009-1336-y.Deichmann, N., T. Kraft, and K. F. Evans (2014), Identification of faults activatedduring the stimulation of the Basel geothermal project from cluster analysis andfocal mechanisms of the larger magnitude events, Geothermics , , 84–97, doi:10.1016/j.geothermics.2014.04.001.Denolle, M. A., E. M. Dunham, G. A. Prieto, and G. C. Beroza (2013), Ground motionprediction of realistic earthquake sources using the ambient seismic field, Journal ofGeophysical Research: Solid Earth , (5), 2102–2118, doi:10.1029/2012JB009603. D R A F T D R A F T - 36
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
Denolle, M. A., E. M. Dunham, G. A. Prieto, and G. C. Beroza (2014), Strong groundmotion prediction using virtual earthquakes,
Science , (January), 399–404, doi:10.1126/science.1245678.Douglas, J., B. Edwards, V. Convertito, N. Sharma, A. Tramelli, D. Kraaijpoel, B. M.Cabrera, N. Maercklin, and C. Troise (2013), Predicting ground motion from inducedearthquakes in geothermal areas, Bulletin of the Seismological Society of America , (3), 1875–1897, doi:10.1785/0120120197.Dyer, B. (2005), Soultz GPK4 stimulation September 2004 to April 2005. Seismic moni-toring report, Semore Seismic Report, Tech. rep.
Dyer, B. C., U. Schanz, T. Spillmann, F. Ladner, and M. O. Haering (2010), Applicationof microseismic multiplet analysis to the Basel geothermal reservoir stimulation events,
Geophysical Prospecting , (5), 791–807, doi:10.1111/j.1365-2478.2010.00902.x.Eberhard, D. A. J., J. D. Zechar, and S. Wiemer (2012), A prospective earthquake forecastexperiment in the western Pacific, Geophysical Journal International , (3), 1579–1592, doi:10.1111/j.1365-246X.2012.05548.x.Edwards, B., T. Kraft, C. Cauzzi, P. Kastli, and S. Wiemer (2015), Seismic monitoring andanalysis of deep geothermal projects in St Gallen and Basel, Switzerland, GeophysicalJournal International , (2), 1020–1037, doi:10.1093/gji/ggv059.Ellsworth, W. L. (2013), Injection-induced earthquakes, Science , , 1225,942–1 –1225,942–7, doi:10.1126/science.1225942.Evans, K. F., A. Zappone, T. Kraft, N. Deichmann, and F. Moia (2012), A survey ofthe induced seismic responses to fluid injection in geothermal and CO2 reservoirs inEurope, Geothermics , , 30–54, doi:10.1016/j.geothermics.2011.08.002. D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 37
Frohlich, C., W. L. Ellsworth, W. A. Brown, M. Brunt, J. Luetgert, T. Macdonald, andS. Walter (2014), The 17 May 2012 M4.8 earthquake near Timpson, East Texas: Anevent possibly triggered by fluid injection,
Journal of Geophysical Research , , 581–593, doi:10.1002/2013JB010755.Gaucher, E., M. Schoenball, O. Heidbach, A. Zang, P. A. Fokker, J.-D. van Wees,and T. Kohl (2015), Induced seismicity in geothermal reservoirs: A review of fore-casting approaches, Renewable and Sustainable Energy Reviews , , 1473–1490, doi:10.1016/j.rser.2015.08.026.Genter, A., K. Evans, N. Cuenot, D. Fritsch, and B. Sanjuan (2010), Contribution ofthe exploration of deep crystalline fractured reservoir of Soultz to the knowledge ofenhanced geothermal systems (EGS), Comptes Rendus Geoscience , (7-8), 502–516,doi:10.1016/j.crte.2010.01.006.Genter, A., N. Cuenot, X. Goerke, B. Melchert, B. Sanjuan, and J. Scheiber (2012),Status of the Soultz geothermal project during exploitation between 2010 and 2012, in Thirty-Seventh Workshop on Geothermal Reservoir Engineering , pp. 1–12.G´erard, A., A. Genter, T. Kohl, P. Lutz, P. Rose, and F. Rummel (2006), The deepEGS (Enhanced Geothermal System) project at Soultz-sous-Forˆets (Alsace, France),
Geothermics , (5-6), 473–483, doi:10.1016/j.geothermics.2006.12.001.Gerstenberger, M. C., and D. A. Rhoades (2010), New Zealand earthquake forecast testingcentre, Pure and Applied Geophysics , (8-9), 877–892, doi:10.1007/s00024-010-0082-4.Giardini, D. (2014), The Swiss Competence Center for Energy Research: Supply of Elec-tricity, in SCCER-SoE Annual Conference Zurich , pp. 1–15.
D R A F T D R A F T - 38
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
Gischig, V. S. (2015), Rupture propagation behavior and the largest possible earthquakeinduced by fluid injection into deep reservoirs,
Geophysical Research Letters , , 7420–7428, doi:10.1002/2015GL065072.Gischig, V. S., and S. Wiemer (2013), A stochastic model for induced seismicity basedon non-linear pressure diffusion and irreversible permeability enhancement, GeophysicalJournal International ,
194 (2) , 1229 – 1249, doi:10.1093/gji/ggt164.Goertz-Allmann, B. P., and S. Wiemer (2013), Geomechanical modeling of induced seis-micity source parameters and implications for seismic hazard assessment,
Geophysics , (1), KS25–KS39, doi:10.1190/geo2012-0102.1.Goodwin, P. (2010), The Holt-Winters approach to exponential smoothing: 50 years oldand going strong, Foresight , , 30–33.Gr¨unthal, G. (1998), European Macroseismic Scale 1998 (EMS-98) , vol. 15, 1–99 pp.,Cahiers du Centre Europ´een de G´eodynamique et de S´eismologie, Luxembourg.Gupta, H. K. (2002), A review of recent studies of triggered earthquakes by artificialwater reservoirs with special emphasis on earthquakes in Koyna, India,
Earth-ScienceReviews , (3-4), 279–310.Gutenberg, B., and C. Richter (1944), Frequency of earthquakes in California, Bulletin ofthe Seismological Society of America , , 185 – 188.Hainzl, S., and Y. Ogata (2005), Detecting fluid signals in seismicity data through sta-tistical earthquake modeling, Journal of Geophysical Research , (B5), B05S07, doi:10.1029/2004JB003247.H¨aring, M. O., U. Schanz, F. Ladner, and B. C. Dyer (2008), Characterisationof the Basel 1 enhanced geothermal system, Geothermics , (5), 469–495, doi: D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 39 > = 5 earthquakes in California, SeismologicalResearch Letters , (1), 78–86, doi:10.1785/gssrl.78.1.78.Hirschberg, S., S. Wiemer, and P. Burgherr (2015), Energy from the Earth Energy fromthe Earth Deep Geothermal as a Resource , 526 pp.Horton, S. (2012), Disposal of hydrofracking waste fluid by injection into subsurfaceaquifers triggers earthquake swarm in Central Arkansas with potential for damagingearthquake,
Seismological Research Letters , (2), 250–260, doi:10.1785/gssrl.83.2.250.Kagan, Y. Y. (2010), Statistical distributions of earthquake numbers: consequenceof branching process, Geophysical Journal International , (3), 1313–1328, doi:10.1111/j.1365-246X.2009.04487.x.Kagan, Y. Y., and D. D. Jackson (2000), Probabilistic forecasting of earthquakes, Geo-physical Journal International , (2), 438–453, doi:10.1046/j.1365-246X.2000.01267.x.Karvounis, D. C., and S. Wiemer (2015), Decision making software for forecasting inducedseismicity and thermal energy revenues in enhanced geothermal systems, in ProceedingsWorld Geothermal Congress 2015 , April, pp. 1–10.Keranen, K. M., H. M. Savage, G. A. Abers, E. S. Cochran, K. M. Keranen, H. M. Savage,G. A. Abers, and E. S. Cochran (2013), Potentially induced earthquakes in Oklahoma,USA: Links between wastewater injection and the 2011 M w 5 . 7 earthquake sequence,
Geology , (March), 1–5, doi:10.1130/G34045.1.Kim, W.-Y. (2013), Induced seismicity associated with fluid injection into a deepwell in Youngstown, Ohio,
Journal of Geophysical Research , , 3506–3518, doi: D R A F T D R A F T - 40
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
Thirty-NinthWorkshop on Geothermal Reservoir Engineering, Stanford University, Stanford, Cali-fornia , pp. 1–9.Kohl, T., and T. M´egel (2007), Predictive modeling of reservoir response to hydraulicstimulations at the European EGS site Soultz-sous-Forˆets,
International Journal of RockMechanics and Mining Sciences , (8), 1118–1131, doi:10.1016/j.ijrmms.2007.07.022.Kraft, T., S. Wiemer, N. Deichmann, T. Diehl, B. Edwards, A. Guilhem, F. Haslinger,E. Kiraly, E. Kissling, A. Mignan, K. Plenkers, D. Roten, S. Seif, and J. Woessner(2013), The ML 3.5 induced earthquake sequence at Sankt Gallen, Switzerland, in Abstract S31F-03, presented at 2013 Fall Meeting, AGU, San Francisco, CA, 9-13 Dec.
Langenbruch, C., and S. A. Shapiro (2010), Decay rate of fluid-induced seismicity aftertermination of reservoir stimulations,
Geophysics , (6), 1–10.Lombardi, A. M., and W. Marzocchi (2010), The assumption of poisson seismic-rate vari-ability in CSEP/RELM experiments, Bulletin of the Seismological Society of America , (5 A), 2293–2300, doi:10.1785/0120100012.Majer, E. L., R. Baria, M. Stark, S. Oates, J. Bommer, B. Smith, and H. Asanuma (2007),Induced seismicity associated with Enhanced Geothermal Systems, Geothermics , (3),185–222, doi:10.1016/j.geothermics.2007.03.003.McClure, M. W. (2012), Modeling and characterization of hydraulic stimulation and in-duced seismicity in geothermal and shale gas reservoirs, Phd, Stanford University. D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 41
McClure, M. W. (2015), Generation of large postinjection-induced seismic events by back-flow from dead-end faults and fractures,
Geophysical Research Letters , , 6647–6654,doi:10.1002/2015GL065028.McClure, M. W., and R. N. Horne (2012), Investigation of injection-induced seismicityusing a coupled fluid flow and rate / state friction model, Geophysics , (6), WC181 –WC198, doi:10.1190/GEO2011-0064.1.McGarr, A. (1976), Seismic moment and volume changes, Journal of Geophysical Re-search , (8), 1478 – 1494.Mena, B., S. Wiemer, and C. Bachmann (2013), Building robust models to forecast theinduced seismicity related to geothermal reservoir enhancement, Bulletin of the Seis-mological Society of America , (1), 383–393, doi:10.1785/0120120102.Mignan, A. (2015), Static behaviour of induced seismicity, Nonlinear Processes in Geo-physics , , 1–16, doi:10.5194/npgd-22-1-2015.Mignan, A., C. Jiang, J. D. Zechar, S. Wiemer, Z. Wu, and Z. Huang (2013), Completenessof the Mainland China earthquake catalog and implications for the setup of the Chinaearthquake forecast testing center, Bulletin of the Seismological Society of America , (2 A), 845–859, doi:10.1785/0120120052.Nanjo, K. Z., H. Tsuruoka, N. Hirata, and T. H. Jordan (2011), Overview of the firstearthquake forecast testing experiment in Japan, Earth, Planets and Space , (3), 159–169, doi:10.5047/eps.2010.10.003.Obermann, A., T. Kraft, E. Larose, and S. Wiemer (2015), Potential of ambient seismicnoise techniques to monitor the St . Gallen geothermal site (Switzerland), Journal ofGeophysical Research : Solid Earth , , 1–16, doi:10.1002/2014JB011817. D R A F T D R A F T - 42
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
Olivella, S., J. Carrera, a. Gens, and E. E. Alonso (1994), Nonisothermal multiphase flowof brine and gas through saline media,
Transport in Porous Media , (3), 271–293,doi:10.1007/BF00613282.Orlecka-Sikora, B. (2010), The role of static stress transfer in mining induced seis-mic events occurrence, a case study of the Rudna mine in the Legnica-Glogow Cop-per District in Poland, Geophysical Journal International , (2), 1087–1095, doi:10.1111/j.1365-246X.2010.04672.x.Reasenberg, P. A., and L. M. Jones (1989), Earthquake Hazard After a Mainshock inCalifornia, Science , , 1173 – 1176.Rhoades, D. A., D. Schorlemmer, M. C. Gerstenberger, A. Christophersen, J. D. Zechar,and M. Imoto (2011), Efficient testing of earthquake forecasting models, Acta Geophys-ica , (4), 728–747, doi:10.2478/s11600-011-0013-5.Rinaldi, A. P., V. Vilarrasa, J. Rutqvist, and F. Cappa (2015), Fault reactivation during CO sequestration: Effects of well orientation on seismicity and leakage, GreenhouseGases: Science and Technology , (5), 645–656, doi:10.1002/ghg.1511.Ruckstuhl, A. (2014), Robust Fitting of Parametric Models Based on M-Estimation, Lec-ture notes , https://stat.ethz.ch/wbl/wbl4/WBL4 robstat14E.pdfRutqvist, J., and O. Stephansson (2003), The role of hydromechanical coupling in frac-tured rock engineering,
Hydrogeology Journal , , 7 – 40.Schoenball, M., C. Baujard, T. Kohl, and L. Dorbath (2012), The role of triggering bystatic stress transfer during geothermal reservoir stimulation, Journal of GeophysicalResearch: Solid Earth , (B9), B09,307, doi:10.1029/2012jb009304. D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 43
Schorlemmer, D., A. Christophersen, A. Rovida, F. Mele, M. Stucchi, and W. Marzocchi(2010), Setting up an earthquake forecast experiment in Italy,
Annals of Geophysics , (3), 1–9, doi:10.4401/ag-4844.Secanell, R., D. Carbon, F. Dunand, and C. Martin (2009), AP5000 Report - SeismicHazard and Risk assessments during three reference time periods (normal, stimulationand circulation), Tech. Rep. October 2009 .Seeber, L., J. G. Armbruster, and W.-Y. Kim (2004), A Fluid-injection-triggered earth-quake sequence in Ashtabula, Ohio: implications for seismogenesis in stable continentalregions,
Bulletin of the Seismological Society of America , (1), 76–87.Segall, P. (1989), Earthquakes triggered by fluid extraction, Geology , , 942 – 946, doi:10.1130/0091-7613(1989)017 < > Journal of Geophysical Research: Solid Earth , , 1–22, doi:10.1002/2015JB012060.Received.Shapiro, S. A., C. Dinske, C. Langenbruch, and F. Wenzel (2010), Seismogenic index andmagnitude probability of earthquakes induced during reservoir fluid stimulations, TheLeading Edge - Special Section: Microseismic , March 2010 , 304–309.Steacy, S., J. Gomberg, and M. Cocco (2005), Introduction to special section: Stress trans-fer, earthquake triggering, and time-dependent seismic hazard,
Journal of GeophysicalResearch , (B5), B05S01, doi:10.1029/2005JB003692.Student (1908), The Probable Error of a Mean, Biometrika , (1), 1 – 25.Taroni, M., J. D. Zechar, and W. Marzocchi (2013), Assessing annual globalM6+ seismicity forecasts, Geophysical Journal International , (1), 422–431, doi: D R A F T D R A F T - 44
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
Chemometrics and Intelligent Laboratory Systems , (2), 127–136, doi:10.1016/j.chemolab.2004.06.003.Wang, X., and A. Ghassemi (2012), A 3D thermal-poroelastic model for geothermal reser-voir stimulation, in , pp. 1–11.Weingarten, M., S. Ge, J. W. Godt, B. A. Bekins, and J. L. Rubinstein (2015), High-rateinjection is associated with the increase in U.S. mid-continent seismicity, Science , ,1336–1340, doi:10.1126/science.aab1345.Zechar, J. D., and T. H. Jordan (2010), Simple smoothed seismicity earthquake forecastsfor Italy, Annals of Geophysics , (3), 99–105, doi:10.4401/ag-4845.Zechar, J. D., D. Schorlemmer, M. Liukis, J. Yu, F. Euchner, P. J. Maechling, and T. H.Jordan (2010a), The Collaboratory for the Study of Earthquake Predictability per-spective on computational earthquake science, Concurrency Computation Practice andExperience , (12), 1836–1847, doi:10.1002/cpe.1519.Zechar, J. D., M. C. Gerstenberger, and D. A. Rhoades (2010b), Likelihood-based testsfor evaluating space-rate-magnitude earthquake forecasts, Bulletin of the SeismologicalSociety of America , (3), 1184–1195, doi:10.1785/0120090192.Zechar, J. D., D. Schorlemmer, M. J. Werner, M. C. Gerstenberger, D. A. Rhoades,and T. H. Jordan (2013), Regional Earthquake Likelihood Models I: First-orderresults, Bulletin of the Seismological Society of America , (2A), 787–798, doi:10.1785/0120120186. D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 45
Zhuang, J., D. Harte, M. Werner, S. Hainzl, and S. Zhou (2012), Basic models of seismic-ity: temporal models,
Community Online Resource for Statistical Seismicity Analysis ,doi:10.5078/corssa-79905851.Zoback, M. L. (1992), First- and second-order patterns of stress in the lithosphere:the world stress map project,
Journal of Geophysical Research , (B8), 11,703, doi:10.1029/92JB00132. Figure 1.
Flowchart of the Adaptive Traffic Light System. GMM stands for Ground MotionModels, w denotes weighting, PSHA means Probabilistic Seismic Hazard Assessment. The dottedgray line delineates the scope of this paper.
Figure 2.
Seismicity of the Basel 2006 ( a. ) and Soultz-sous-Forˆets 2004 ( b. ) geothermalproject in NS cross sections. Wells are represented by black lines. Dotted light gray gridsindicate voxels of 200 m × m × m for testing. Colors denote moment magnitudes of theevents, note the different scales. Map inset shows the location of the geothermal sites. D R A F T D R A F T - 46
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
Figure 3. A.
Explanation of the spatial component of the SaSS model. a. General equation ofsmoothed seismicity. b. Seismicity of a learning period, colors denote temporal weights indicatedin the inset. c.
3D Gaussian kernel represented in 2D with gray dashed and solid black lines. d. Vertical cross section of smoothed seismicity, colors denote the spatial probability densityfunction. B. Explanation of the HySei model. Black circles represent simulated seismicity on thefault plane. Red dots indicate observed seismicity of the learning period. Minimum, maximum,and intermediate principal axes are a z , a x , and a y respectively. Solid black lines indicate thelength of the principal axes. Black ellipse shows the 95% of the seismicity cloud. Red curve showsthe empirical event distribution along the minimum principal axis (i.e., off-plane direction). C. Explanation of time periods used for model calibration and forecasts.
Figure 4.
Cross sections of the SaSS (upper panels) and HySei (lower panels) forecasts forthe period containing the largest event in the sequence a few hours after the shut-in (2006-12-0816:48, m w . D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 47
Figure 5.
Forecast of the number of events with different learning periods. a. Evolutionof injection flow rate in l/s (black line) and wellhead pressure in MPa (orange line) during theinvestigated time period of the Basel 2006 project. Dotted gray line indicates the shut-in. Blue,red, green and purple horizontal lines correspond to the learning periods starting at the beginingof the injection, ending after 1 .
25 days, 3 .
25 days, 5 .
25 days, and 9 .
50 days, respectively. b. Number of events in 6-hour time bins in function of time for SaSS model. Black dots represent theobserved seismicity rate. Blue, red, green and purple vertical lines correspond to the previouslymentioned learning periods. Shaded areas indicate the corresponding 72-hour forecasts with95% Poissonian confidence intervals. Dotted gray line indicates the shut-in. c. Same as b. forHySei model. d. Same as a. for the Soultz-sous-Forˆets 2004 project with different moments offorecasts: blue, red, green, and purple lines correspond to the learning periods of day 1 .
75, 3 . .
5, respectively. e. Same as b. for the Soultz-sous-Forˆets 2004 project. f. Same as c. for the Soultz-sous-Forˆets 2004 project. D R A F T D R A F T - 48
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
Figure 6.
Number-test. Horizontal axes denote the length of learning periods used forforecasting, vertical axes denote individual 6-hour forecast time windows. a. Observed seismicityrate of the Basel 2006 project. Color corresponds to number of observed events in 6-hour timebins. Dotted red line indicates the shut-in moment. Blue, red, green and purple rectanglesindicate forecasts (vertical direction) with learning period of 1 .
25 days, 3 .
25 days, 5 .
25 days and9 .
50 days, respectively. These forecasts are explicitly plotted on the previous figure. b. Differencebetween number of forecast events by SaSS model and the number of observed events. Blue andred show moments when models overestimate and underestimate the observed seismicity rate,respectively. Yellow indicates similar number of forecast events as observed earthquakes. Grayupward-pointing triangles and solid downward-pointing denote moments when the number offorecast events (with Poissonian error bars) are significantly higher or lower than the numberof observed seismicity rate, respectively. Dotted black line indicates the shut-in moment. c. Same as b. for HySei model. d. Same as a. for the Soultz-sous-Forˆets 2004 project. Note thedifferent color scale. e. Same as b. for the Soultz-sous-Forˆets 2004 project. f. Same as c. forthe Soultz-sous-Forˆets 2004 project. D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 49
Figure 7.
Snapshots of magnitude frequency distribution of observed and forecast earthquakes(forecast is normalized so that total number of forecast earthquakes is equal to the number ofobserved events). Solid squares denote observed seismicity, colors refer to the same moments ason the previous two figures. Orange dots show SaSS forecast rate, blue dots indicate the HySeiforecast rate. Corresponding transparent shaded areas indicate the 95% confidence interval ofthe forecasts. Greenish shaded area indicates the overlapping of the two confidence intervals. a. Seismicity of the 3-day forecast period after 1 .
25 day of learning period. b. Same as a. after 3 . c. Same as b. after 5 .
25 days of learning period. d. Same as a. after 9 . e. Same as a. after 1 .
75 days of learning period in Soultz-sous-Forˆets2004. Note that last magnitude bin contains all magnitudes higher than 1 . f. Same as e. after3 .
50 days of learning period. g. Same as e. after 5 .
00 days of learning period. h. Same as e. after 6 .
50 days of learning period.
Figure 8.
Space-test. Colorbar indicates joint log-likelihood values; yellow indicates betterforecasts than red. Crosses represent moments when the model does not pass the Space-test.Gray squares denote moments when no earthquake occurred. Gray dotted line marks the shut-inmoment. a. Space-test of SaSS for Basel 2006. b. Space-test of HySei for Basel 2006. c. Sameas a. for Soultz-sous-Forˆets 2004. d. Same as b. for Soultz-sous-Forˆets 2004. D R A F T D R A F T - 50
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
Figure 9.
Comparisons of model components based on log-likelihood. Green indicates mo-ments when SaSS is superior to HySei, red shows moments when HySei performs better thanSaSS. White denotes moments when both models perform similarly. Dotted black line indicatesthe shut-in moment. Note that scales are different for each component. a. Comparison of num-ber components for Basel 2006. b. Comparison of magnitude components for Basel 2006. c. Comparison of spatial components for Basel 2006. d. Same as a. for Soultz-sous-Forˆets 2004. e. Same as b. for Soultz-sous-Forˆets 2004. f. Same as c. for Soultz-sous-Forˆets 2004. Figure 10.
Cumulative joint LL values of the entire model summed for the indicated forecastperiods and divided by the total number of observed events in the given forecast period. SaSSis denoted by solid orange dots, HySei is shown by blue circles. Dotted black lines show theshut-in moment. All values are plotted at the end of the corresponding learning period. Regime A indicates the period when SaSS and HySei perform similarly, regime B indicates the periodwhen SaSS performs better than HySei, regime C indicates the period when HySei performsbetter than SaSS. a. Cumulative joint LL/Eqk for 6-hour forecast periods of the Basel 2006experiment. b. Same as a. for 24-hour forecast periods. c. Same as a. for 48-hour forecastperiods. d. Same as a. for 72-hour forecast periods. e. Same as a. for Soultz-sous-Forˆets 2004. f. Same as b. for Soultz-sous-Forˆets 2004. g. Same as c. for Soultz-sous-Forˆets 2004. h. Sameas d. for Soultz-sous-Forˆets 2004. D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 51
Figure 11.
Comparison of different methods to evaluate average information gain (HySeiis the reference model). Solid black line shows the empirical probability distribution of theinformation gain values, dotted gray line indicates the normal distribution, which expected valueand standard deviation is calculated from the information gain population. Dashed red, dashedgreen, solid orange, and solid brown lines correspond to classical mean, robust mean, bootstrapmean, and bootstrap median, respectively. Insets show the estimated average values with theiruncertainties. Reddish background denotes the area, where SaSS is better, yellowish backgrounddenotes the area, where HySei is better. The number of investigated earthquakes are shown inthe top left corner of the graph. a. Information gain in the case of the Basel 2006 experiment. b. Same as a. for Soultz-sous-Forˆets 2004. D R A F T D R A F T - 52
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS ... ww Figure 1.
D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 53
Magnitude Magnitude a. Basel 2006
NS coordinates [m] -2400 -2000 -1600 -1200-5400-5000-4600-4200
NS cross sectionNS cross section Distance from the GPK1 wellhead in NS [m] D e p t h [ m ] D e p t h [ m ] b. Soultz-sous-Forêts 2004
Basel (CH)Soultz-sous-Forêts (FR)
Figure 2.
D R A F T D R A F T - 54
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS Σ All observed events
Past seismicity with temporal weighting 3D Gaussian kernel Smoothed seismicity x [m]-500 500-500700 d. z [ m ] z [ m ] x [m]-500700 -300 700 Σ b. c. All observed events
Distance from (x eqk , z eqk ) [m] W e i g h t i n s p a c e -dx;-dz01 dx;dz00 1Weight 0 4x10 -2 PDF a. W e i g h t Days0 8
A.C.B.
SaSS modelHySei modelTime periods
Long es t p r i n c i p l e ax i s [ m ] -600 -400 -200 0 200 400 600 -600 6004002000-200-400 -300-100100300 a a z a x Empirical distributionof observed events along a z I n t e r m e d i a t e p r i n c i p l e a x i s [ m ] Events simulated by HySei S ho r t es t p r i n c i p l e ax i s [ m ] a x a y Figure 3.
D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 55 P D F F i g u r e . D R A F T D R A F T - 56
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
LP1 LP2 LP3 LP4 O b se r ve d se i s m i c i t y L P = . d L P = . d L P = . d L P = . d S hu t- i n LP1 LP2 LP3 LP4 B ase l L P L P L P L P a . b . c . S a SS m od e l H y d r a u li cs H y S e i m od e l E l a p se d d ays f r o m t h e s t a r t o f t h e s t i m u l a t i on Number of events / 6hNumber of events / 6hInjected flow rate [l/s] Wellhead pressure [MPa]
LP1 LP2 LP3 LP4 L P d . e . f . S a SS m od e l H y d r a u li cs H y S e i m od e l E l a p se d d ays f r o m t h e s t a r t o f t h e s t i m u l a t i on Number of events / 6hNumber of events / 6h S ou l t z - s ou s - Fo r ê t s LP1 LP2 LP3 LP4 O b se r ve d se i s m i c i t y L P = . d L P = . d L P = . d L P = . d S hu t- i n L P L P L P Injected flow rate [l/s] Wellhead pressure [MPa] F i g u r e . D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 57 S a SS m od e l H y S e i m od e l L e ng t h o f l ea r n i ng p e r i od [ d ays ] Forecast time window [days] L e ng t h o f l ea r n i ng p e r i od [ d ays ] N u m b e r o f ob se r ve d ea r t hqu akes a . b . c . N u m b e r o f ob se r ve d ea r t hqu akes d . S a SS m od e l e . o v e r p r ed i c t unde r p r ed i c t pe r f e c tf o r e c a s t -
50 0 50 D i ff H y S e i m od e l f . -
50 0 50 o v e r p r ed i c t unde r p r ed i c t pe r f e c tf o r e c a s t D i ff B ase l S ou l t z - s ou s - Fo r ê t s F i g u r e . D R A F T D R A F T - 58
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
Magnitudes N u m b e r o f eve n t s Forecast from day 1.75 for 72hForecast from day 5.00 for 72h Forecast from day 6.50 for 72h
Magnitudes N u m b e r o f eve n t s Basel 2006Soultz-sous-Forêts 2004 e. f.g. h. N u m b e r o f eve n t s Forecast from day 1.25 for 72ha. b. Forecast from day 9.50 for 72hForecast from day 5.25 for 72h
Magnitudes d. Forecast from day 3.25 for 72h
Magnitudes N u m b e r o f eve n t s c. Forecast from day 3.50 for 72h Observed RateSaSS Forecast RateHySei Forecast RateObserved RateSaSS Forecast RateHySei Forecast Rate
Figure 7.
D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 59 H y S e i m od e l H y S e i m od e l L e ng t h o f l ea r n i ng p e r i od [ d ays ] L e ng t h o f l ea r n i ng p e r i od [ d ays ] Forecast time windows [days] B ase l S ou l t z - s ou s - Fo r ê t s - - -
15 a . b . c . d . S a SS m od e l S a SS m od e l LL F i g u r e . D R A F T D R A F T - 60
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS LL S i S - LL H y S e i L e ng t h o f l ea r n i ng p e r i od [ d ays ] L e ng t h o f l ea r n i ng p e r i od [ d ays ] Forecast time windows [days] SpatialMagnitude Number B ase l S ou l t z - s ou s - Fo r ê t s H y S e i i s b e tt e r S a SS i s b e tt e r - - - . . H y S e i i s b e tt e r S a SS i s b e tt e r a . b . c . d . e . f . H y S e i i s b e tt e r S a SS i s b e tt e r - F i g u r e . D R A F T D R A F T
IRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
X - 61
Soultz-sous-Forêts 2004 -12-10-8-6-4-20 -12-10-8-6-4-20 -12-10-8-6-4-20 -15-10-50 h f o r ecas t h f o r ecas t h f o r ecas t h f o r ecas t j LL / E q k -30-25-20-15-10-50 -25-20-15-10-50 -20-15-10-50 HySeiSaSS A B C j LL / E q k j LL / E q k j LL / E q k Length of learning period [days]Basel 2006 a. Length of learning period [days]
B C b.c.d. e.f.g.h. HySeiSaSS
Figure 10.
D R A F T D R A F T - 62
KIRALY-PROAG ET AL.: VALIDATING INDUCED SEISMICITY FORECAST MODELS
Basel 2006
Information gain
Soultz-sous-Forêts 2004 -4 -2 0 2 4 6 8 10 1200.10.20.30.40.50.60.7 P r ob a b ili t y d e n s i t y Information gain P r ob a b ili t y d e n s i t y -4 -2 0 2 4 6 8 10 1200.10.20.30.40.50.60.70.80.9 Information gain: SaSS / HySei IG SaSS is better HySei is betterClassical meanBootstrap mean Robust meanBootstrap median SaSS is better n eqk = 1695n eqk = 900 Average information gain
Classical meanBootstrap mean Robust meanBootstrap median
Average information gaina. b.