A stochastic model for induced seismicity based on non-linear pressure diffusion and irreversible permeability enhancement

During deep reservoir engineering projects, in which permeability is enhanced by high-pressure ﬂuid injection, seismicity is invariably induced, posing nuisance to the local population and a potential hazard for structures. Hazard and risk assessment tools that can operate in real-time during reservoir stimulation depend on the ability to efﬁciently model induced seismicity. We here propose a novel modelling approach based on a combination of physical considerations and stochastic elements. It can model a large number of synthetic event catalogues, and at the same time is constrained by observations of hydraulic behaviour in the injection well. We model ﬂuid ﬂow using non-linear pressure diffusion equations, in which permeability increases irreversibly above a prescribed pressure threshold. The transient pressure ﬁeld is used to trigger events at so-called ‘seed points’ that are distributed randomly in space and represent potential earthquake hypocentres. We assign to each seed point a differential stress based on the mean estimates of the in situ stress ﬁeld and add a normal distributed random value. Assuming a fault orientation with respect to the stress ﬁeld and a Mohr–Coulomb failure criterion, we evaluate at each time step, if a seed point is triggered through a pressure increase. A negative proportional relationship between differential stress and b values is further assumed as observed from tectonic earthquakes and in laboratory experiments. As soon as an event is triggered, we draw a random magnitude from a power-law distribution with a b value corresponding to the differential stress at the triggered seed point. We thus obtain time-dependent catalogues of seismic events including magnitude. The strategy of modelling ﬂow and seismicity in a decoupled manner ensures efﬁciency and ﬂexibility of the model. The model parameters are calibrated using observations from the Basel deep geothermal experiment in 2006. We are able to reproduce the hydraulic behaviour, the space-time evolution of the seismicity and its frequency–magnitude distribution. A large number of simulations of the calibrated model are then used to capture the variability of the process, an important input to compute probabilistic seismic hazard. We also use the calibrated model to explore alternative injection scenarios by varying injection volume, pressure as well as depth, and show the possible effect of those parameters on seismic hazard.


I N T RO D U C T I O N
The Earth's crust is critically stressed in most places, and it is a well established fact that man-made perturbations to the stress conditions-for example through reservoir stimulation, hydraulic fracturing (e.g. Rutledge et al. 2004), mining (e.g. Simpson 1986;Kwiatek et al. 2010), reservoir impoundment (e.g. Talwani 1997;Gupta 2002;Chen 2009), injection of waste water (e.g. Nicholson & Wesson 1990) or CO 2 storage (e.g. Rutqvist et al. , 2011aNichol et al. 2011)-leads to enhanced seismic activity. Induced seismicity has received increased attention in the past few years, because a number of GeoEnergy projects have been delayed or cancelled because of felt earthquakes and the concern they caused.
Although the magnitude of such induced seismic events is typically smaller than the largest observed natural events in the same location, they are governed by the same physics and are often indistinguishable from natural events (e.g. . Although induced seismicity is related to numerous human activities, it has recently been studied most extensively in the context of the creation of enhanced geothermal systems (EGS). The potential of EGS in contributing to the future energy strategy has long been recognized (Tester 2006). They are a particularly attractive energy system, because they can be principally installed in most locations, regardless of the presence of a favourable hydrothermal system. Induced seismicity during reservoir stimulation is one of the major obstacles for standard use of EGS near urban centres.
Although stimulation with fluids is essential for enhancing permeability and thus enabling productivity of the reservoir, the required shearing and dilating of fractures is associated with seismic energy radiation. Induced seismic events with magnitudes exceeding M3 potentially represent a risk for infrastructure, in extreme cases even for people. However, even smaller but felt events may already have a severe impact on public acceptance. The reaction to nuisance and fear conjured by felt seismic events led to suspension of the planned EGS at Basel (Häring et al. 2008), and the formation of citizen groups in other instances. However, induced seismicity is both a blessing and a curse: it also gives important insights into size of the created reservoir. Thus, EGS experiments are typically monitored with down-hole instrumentation, and the hundreds to thousands of events recorded are a rich data set to calibrate models of hydromechanical interaction. Majer et al. (2012) define steps how to address induced seismicity, and express the need of quantifying hazard and risk related to it. However, reliable and validated methods to do so do not exist to date. Managing hazard and risk of induced seismicity during reservoir stimulation requires recognizing potentially hazardous levels of seismicity early, ideally in a risk study conducted even before the expensive drilling phase has started. To date, however, there are no reliable predictors of the response of the underground at a given location to injection. Therefore, meaningful hazard and risk assessment of such projects must be mainly conducted during the ongoing reservoir creation, so that appropriate actions can be taken before harmful events occur. The method used so far in EGS projects is the so-called traffic light system, in which different sets of actions are defined during project planning based on thresholds of peak-ground velocity, magnitude and public perception [e.g. used in Berlìn, El Salvador (Bommer et al. 2006) and in Basel (Häring et al. 2008)]. The largest deficit of the method is that actions are taken only if predefined thresholds are exceeded, that is it is only reactive. If instead seismicity rates depending on time and magnitude, including its inherent uncertainties, can be modelled up to a certain time during injection, a probabilistic forecast of seismicity, hazard or risk occurring during the subsequent hours can be computed. Decisions on future injection strategies (e.g. increasing or reducing injection rate, shut-in or bleed-off of the well) can then be taken based on a probabilistic threshold of a certain event to occur (e.g. magnitude, ground motion, monetary loss, etc.), already before a harmful event happens. Such a hazard management scheme based on probabilistic seismic hazard analysis (PSHA) relies on the ability to model time-dependent seismic hazard. For this both validated statistical methods and numerical reservoir simulators are urgently needed. The conceptual idea and current state of development of approaches in both model classes are described in the following.

Statistical forecast models
Statistical methods to forecast seismic hazard are described by Bachmann et al. (2011) and Mena et al. (2013), and have been tested in a pseudo-prospective manner for the EGS stimulation at Basel in 2006. Because of their robustness and efficiency, they are most useful for real-time traffic light applications. Mena et al. (2013) show that both an adaptation of a widely used earthquake clustering model, the Epidemic Type Aftershock Model (Ogata 1992;Bachmann et al. 2011), as well as a model suggested by Shapiro et al. (2010) can forecast the seismicity during the Basel stimulation quite well in a simulated near-real-time application. The Shapiro model forecasts the rate of induced seismicity during injection experiments as a function of a site-specific parameter (the so-called seismogenic index) and of the injected fluid volume. Mena et al. (2013) also show that a combination of these models, were the relative weights are updated at each time step according to the relative performance of each model, is superior to an individual model in its robustness and forecasting ability. However, while powerful in near-real-time applications, such statistical models have clear limits: They account for the underlying physical processes governing induced seismicity only to a limited degree. The two aforementioned models only consider injection volume as information. All other parameters rely on calibration against observed seismicity data, which has to be done for each individual site. Hence, they have limited forecasting capabilities for longer time periods, alternative injection scenarios, and post-shut-in behaviour.

Physics-based models
Physics-based models can potentially overcome the shortcomings of statistical models that only include a limited degree of physics. As they describe physical processes more comprehensively, they may perform better as forecasts model. In addition, they have the capability of exploring the sensitivity of reservoir performance and induced seismicity to various processes (e.g. stress redistribution, thermal contraction, etc.), site specific conditions (e.g. initial hydraulic properties, in situ stress state, etc.) and design parameters that can be controlled in a project (e.g. injection volume or pressure, reservoir depth and size). However, a drawback of those models is the numerous parameters that are often badly constraint and difficult to calibrate against observations. Such physics-based models rely on numerical methods to simulate (thermo-)hydromechanical processes in geothermal reservoirs. So far they have mostly been used in scenario-type applications, because running them in near-real time during an ongoing stimulation is challenging. Generally, full thermohydromechanical coupling in a fractured medium containing multiphase fluids has to be included in a numerical model to appropriately account for most phenomena associated with fluid-driven seismicity. Another essential but rarely met requirement for their use in a seismic hazard analysis framework is the ability to forecast the magnitude distribution of induced events. We here give a short review of simulators that have been used in the induced seismicity context: The code FRACAS was presented by Cacas et al. (1990) and Bruel (2007), and applied to the EGS system at Soultz-sous-Forêt (Baujard & Bruel 2006). The predefined fractures are assigned a stress state depending on the tectonic stress field and a failure criterion that governs the ability of the fractures to shear. Fracture permeability is updated as a function of shearing-induced dilation. The code HEX-S presented by Kohl & Mégel (2007) similarly accounts for enhancements of fracture permeability through shear dilation. The model does not include stress transfer and cannot compute event magnitudes. A 2-D model was developed by Baisch et al. (2010) to account for both slip-dependent permeability and stress transfer within the modelling plane. It can calculate the magnitude of individual events from the amount of slip and the slipped area. It was applied in the risk study conducted after the Basel injection experiment (Baisch et al. 2009a). It was able to reproduce a number of characteristics of the induced seismicity, such as post-injection seismicity and the occurrence of the largest events after shut-in. McClure & Horne (2011) present an approach which includes slip governed by rate-and-state friction, and present a generic study to explore the role of injection pressure on magnitudes of induced events. They found that larger injection pressure results in larger magnitudes. Rutqvist et al. (2002) and Rutqvist (2011) suggest a combination of the far-developed commercial simulators FLAC3D and TOUGH2. Pressure diffusion and heat transport are solved with TOUGH2. At each time step the pressure and temperature fields are transferred to FLAC3D, which solves the hydrothermomechanical response of the rock mass. The method was applied at The Geysers geothermal site to explore effects of decreasing reservoir temperatures on enhanced seismicity (Rutqvist & Oldenburg 2008). The approach is very powerful in including various processes associated with fluid properties and reservoir mechanics. However, it currently cannot simulate the magnitude of large numbers of earthquakes.
Although the aforementioned physics based approaches are successful in simulating various phenomena associated with reservoir creation and induced seismicity, none of them has ever been used in induced seismicity PSHA. Most existing models allow stochastic variability of parameters (such as fracture orientation or extent, friction parameters, stress parameters, etc.), but they do not systematically present the uncertainty in our knowledge of these critical parameters. They also typically do not forecast meaningful magnitude distributions that extrapolate to the very rarely observed events. These events, despite being rare, dominate the hazard and risk at lower probability levels (e.g. Mena et al. 2013). Bachmann et al. (2012) and Goertz-Allman & Wiemer (2013) introduced a so-called hybrid model to be used in PSHA, which strives to combine the advantages of statistical and physical models. They use a simple linear flow model and a 'stochastic seed model' build from basic geomechanical considerations similar as suggested by Rothert & Shapiro (2003). In addition to the approach by Rothert & Shapiro (2003), it includes a method for producing magnitudes. It can thus be calibrated against observed seismicity, for example by adjusting to the density of faults in the stimulated volume. The model is able to explain a wide range of observations: The overall earthquake size distribution (or b value), the observed spatial distribution of b values, the observed stress drop as a function of distance (Goertz-Allmann et al. 2011) and the fact that largest events may often be observed shortly after shut-in. However, the model presented by Goertz-Allmann & Wiemer (2013) is clearly limited by the overly simplistic pressure diffusion model.

Towards hybrid models
In our assessment, the next generation of modelling tools needed for induced seismicity in general and for risk management of EGS systems in particular need to have the following key characteristics: (1) They need to be applicable in both near-real-time applications as one element of advanced traffic light systems, as well as in scenario modelling the design of reservoir stimulation and sensitivity studies.
(2) They need to take into account the inherent uncertainty in our knowledge of critical model parameters, and in our limited understanding of the physical processes at work. A large number of model realizations (several thousands) have to be computed to cover the entire model uncertainty, which is very demanding in terms of computational costs.
(3) They need to include ability to couple both induced seismicity as well as the evolution of reservoir properties based on observations and a meaningful physical representation of the system. This will eventually lead to the ability to manage the inherent trade-off between permeability creation and seismic hazard.
We here present an extension of the hybrid approach by Goertz-Allmann & Wiemer (2013) that can forecast induced seismicity in near-real-time during future injection experiments, and is able to model scenarios for planning geothermal reservoirs before injection. The model represents the underlying key physical processes in a simplistic manner to ensure computation efficiency. The flow simulator accounts for non-linear pressure diffusion with irreversible permeability enhancement. The seismicity component computes synthetic catalogues of induced events including magnitudes based on a stochastic process. Thus, the flow and the seismicity simulators are fully decoupled. This reduces the number of free parameters and enhances efficiency and flexibility of the model. The approach is designed to reproduce the main observables during reservoir stimulation, namely injection rate and pressure, temporal and spatial evolution of statistical seismicity characteristics.
The paper is organized as follows: The next section describes the non-linear flow model and summarizes the stochastic seed model. In Section 3, we describe the approach to calibrate our model parameters using observed data. Results of the calibration procedure using the Basel data set are presented in Section 4. The calibrated model is used to estimate seismic hazard using ground motion prediction equations (GMPE, Section 5), and to explore different injection scenarios and their effect on seismic hazard (Section 6). In Section 7, we discuss the limitations of the model and the possible measures for limiting seismicity in future stimulation experiments. Section 8 summarizes the main results.

M O D E L D E S C R I P T I O N
The workflow of the model is outlined in Fig. 1. We use a flow model calibrated against pre-stimulation and stimulation hydraulic data to compute the transient pressure distribution. The pressure model is then used to decide whether events are triggered at so-called seed points. From a large number of synthetic event catalogues (here 1000), seismicity rates per magnitude bin are derived and used to compute seismic hazard expressed in terms of an EMS intensity (European Macroseismic Scale, EMS98; Gruenthal 1998). In the following sections, each model component is described in detail.

Pore pressure diffusion model
A minimum requirement for meaningful modelling of flow during stimulation, with the simplification of neglecting hydromechanical coupling, is to account for the dramatic and mostly irreversible permeability increase as shear dilation occurs. Observations in the laboratory (Mitchell & Faulkner 2008), during natural earthquake sequences (Miller et al. 2004;Cappa et al. 2009) and during reservoir stimulation experiments (e.g. at Soultz-sous-Fôret, Evans 2005) reveal that permeability increases by two to three orders of magnitude during shearing. Using variable permeability in pressure diffusion models has a large impact on the pressure distribution. Murphy et al. (2004) show that flow in fractures with pressure-dependent aperture results in a pressure field expanding in a shock-wave manner rather than in a smooth diffusive way. At distances on the order of several 10-100 m from the injection point such non-linear diffusion models exhibit pressures that are several orders of magnitudes higher compared to linear diffusion models. Shapiro & Dinske (2009), Hummel & Müller (2009 and Hummel & Shapiro (2012) used non-linear pressure diffusion models with pressure-dependent diffusivity for seismicity-based characterization of reservoir hydraulic parameters as suggested by Shapiro et al. (2002). Ortiz et al. (2011) use an explicitly discretized fracture model combined with a continuum model (using the finite-element modelling package COMSOL) in which the fracture permeability is pressure dependent. Their model is able to reproduce the prestimulation test and the early times of stimulation pressure curve recorded at Basel. The flow model deployed here also includes pressure dependent permeability changes. Unlike to the aforementioned non-linear pressure diffusion models, permeability increases irreversibly above a certain pressure threshold to more realistically account for permeability evolution during stimulations (e.g. Evans 2005). We use the COMSOL simulation software due to its flexibility to combine partial differential equations.

Governing equations
The pressure diffusion eq. (1) is coupled with an expression for the rate of change of a variable u (2), here denoted 'stimulation factor'. The pressure diffusion equation is where ρ (kg m −3 ) is fluid density, S (Pa −1 ) is the specific storage coefficient, κ (m 2 ) is permeability, η is fluid viscosity (Pa s) and q m is a mass source (kg m −3 s −1 ). Because permeability κ is not constant during stimulation, we replace it by with κ 0 being the initial permeability before the stimulation, which is increased using the stimulation factor u. u is calculated as follows: The scalar functions H xx are smoothed Heavyside functions with the following characteristics: H pt is one if pressure increases and zero otherwise, H u is one if u is below the maximum stimulation factor u t . H p is one if pressure is above the threshold pressure p t . Thus, the stimulation factor u (and therefore permeability) increases only if pressure increases, but stays at its current value if pressure decreases again. This mimics the irreversible manner that permeability changes due to shear dilation, as observed in laboratory tests (e.g. Esaki et al. 1999;Lee & Cho 2002) and in EGS experiments (e.g. in Soultz; Evans 2005). In reality, a small portion of fracture dilation is reversible due the fracture compliance; however, we ignore this reversible part of permeability change for simplicity. Furthermore, permeability only increases above a certain pressure threshold p t , here denoted 'stimulation pressure', and below a limiting stimulation factor u t , hereafter called 'stimulation limit'. The latter accounts for the fact that dilation during shearing does not increase infinitely, but converges towards a maximum level as shearing continues (Lee & Cho 2002). Note that H u and H p are smoothed heavy-side functions with transition zones u t and p t about their threshold values (Fig. 2). On the one hand, this ensures numerical stability; on the other hand, it is more realistic as the rate of permeability increase does not start abruptly but does so rather gradually. In all models shown hereafter, u t is arbitrarily set equal to u t. . (Note that smaller values of u t do not have a large impact on the solutions presented later, it just has to be large enough for numerical stability). p t is an adjustable parameter (later calibrated along with the other parameters), and (p t − p t ) defines the onset pressure for stimulation. The factor C u is a constant that scales the rate at which permeability changes and is here referred to as 'stimulation velocity'.

Model geometry, mesh and boundary conditions
The flow models presented here are 2-D, and exhibit radial symmetry. The assumption is justified to some degree, because the seismicity clouds observed at different stimulation experiments [e.g. Basel (Häring et al. 2008), Cooper Basin, Australia (Baisch et al. 2006(Baisch et al. , 2009b, etc.] often exhibit strongly oblate geometries. Thus, both flow and seismicity is modelled only in the plane of maximum extent of the observed seismicity cloud. It can be interpreted as a layer that includes a fault or a fracture zone with limited width. We acknowledge that the 2-D radially symmetric model is only an equivalent continuum representation of the complex fractured nature of the reservoir. We further take advantage of the symmetries of a radially symmetric model, and model only one quadrant of the plane (Fig. 2d). Thus, two boundaries of the model domain are symmetry boundaries and the other two pressure boundary conditions. Because we only model pressures above formation pressure (i.e. pressure perturbations), we set both initial pressure and pressure boundary conditions to 0 MPa. Gravity effects on pressure diffusion are neglected. A triangular mesh with a side length of 15 m in the far field and about half the effective borehole radius close to injection is chosen. Special attention is given to the mass source boundary condition that is a circular hole in the modelling domain representing the borehole. Wellbore pressure measured during hydraulic tests is strongly affected by conditions of the borehole wall (Cinco-Ley & Samaniego 1981; Kruseman & de Ridder 1994;Economides & Nolte 2000). A zone of low permeability close to the wellbore forces the injection pressure at a given flow rate to increase, whereas a strongly damaged (i.e. fractured) zone with high permeability results in lower pressure. Due to the strong impact of such nearwellbore effects on the measured injection pressure curve, they have to be considered to some degree. Modelling near-wellbore effects by explicitly including them into the geometry (i.e. different zones with strong permeability contrast) is very complex and computationally expensive. However, in case of a heavily fractured zone they can also be accounted for by using an effective borehole radius r eff that is larger than the real one (e.g. Kruseman & de Ridder 1994). Thus, we here minimally account for the negative skin effect by applying fluid injection along the boundary of a circle with effective borehole radius r eff . In our case, the effective radius is larger than the actual borehole radius because it represents a heavily fractured zone where the fluid enters the rock mass. The injected fluid mass is uniformly distributed at the boundary. r eff is estimated as explained in the next section.

Model ambiguity
Modelling the transient pressure field of a reservoir based on pressure and flow data of a single borehole leads inevitably to an underdetermined problem and a resulting ambiguity in the estimated hydraulic parameters. Only the permeability κ can be determined unambiguously. We illustrate this ambiguity with a simple example based on a linear flow model for simplicity (i.e. permeability κ is constant), and using a synthetic injection schedule (Fig. 3a). Permeability was set to κ 0 = 10 −15 m 2 . Fig. 3 illustrates the performance of the model for a number of combinations of effective borehole radii r eff and storage coefficients S. Although the pressure curve extracted at the mass boundary ( Fig. 3b) is the same for all combinations of r eff and S shown in Fig. 3(c), the pressure profiles ( Fig. 3d) are fundamentally different. Pressure perturbations reach farther into the reservoir for small storage coefficients and large effective borehole radii. For resolving S and r eff unambiguously, additional information would be required, ideally from the pressure response at a second, nearby well. In absence of such information, the extent of induced seismicity can provide additional constraints on the reservoir hydraulics in case of reservoir stimulation (Shapiro et al. 2002;Shapiro & Dinske 2009). In Section 4, we use the extent of the seismic cloud to constrain S and r eff for the case of the Basel stimulation.

Performance of the non-linear flow model
In Fig. 4, we choose the synthetic injection schedule defined in Fig. 3(a) to illustrate the performance of our non-linear model. As the model results here only serve as explanation, the stimulation parameters are chosen arbitrarily and are identical to the ones used in Fig. 2, (i.e. p t = 4 MPa, u t = 10, u t = u t and p t = 1 MPa) with the constant C u = 0.005 s −1 . Best-fitting parameters later derived for the Basel case study (Section 4) are given in Table 1.
The stimulation factor u increases as long as pressure increases and stays constant otherwise (Figs 4a and b). Note that the maximum possible value for u is not u t, but almost twice its value. This is a consequence of the transition zone u t , which causes ∂u/∂t to gradually decay above u t and becomes zero if u = u t + u t (see Fig. 2b). Fig. 4(c) shows u plotted against pressure, and illustrates the hysteretic nature of irreversible permeability changes. After stimulation, the medium contains a zone of irreversibly increased permeability (Fig. 4d). In Figs 4(e) and (f) we illustrate the role of the stimulation threshold p t by varying it from 2 to 6 MPa. The pressure above which permeability increases affects the distance that pressure disturbances reach. If permeability increases already at low pressures, the pressure front reaches farther from the injection location. Thus, in the case of non-linear pressure diffusion, not only initial permeability and storage [i.e. diffusivity defined as D 0 = κ 0 /(ηS)], defines the extent of pressure disturbance in the medium, but the threshold above which the medium is stimulated has also a major impact. Note that the non-linear model is able to compute a steep pressure front and pressure magnitudes in the reservoir on the order of several MPa, which is in agreement with triggering pressure magnitudes suggested by Terakawa et al. (2012).

The geomechanical stochastic seed model (GSSM)
The described pressure model is a pre-requisite input for the GSSM, which is based on transient changes of effective stress at locations of potential earthquakes. The GSSM is described in detail by Goertz-Allmann & Wiemer (2013). Below we present a short summary and adoption to our case study, along with the results of a sensitivity analysis. Best-fitting parameter values later found using the Monte Carlo approach (Section 3) for calibration against seismicity data recorded at Basel are given in Table 2.

Model set-up
Potential earthquake locations are represented as points uniformly random distributed over the modelling domain (Fig. 5). The model domain corresponds to four mirrored quadrants of the pressure field model. They are called seed points or seed faults. The number of seed points per unit area or volume in the model is the seed density N s . For each seed point, a maximum and minimum stress magnitude σ 1 and σ 3 are drawn from a normal distribution with the average taken from in situ stress estimates [taken from Häring et al. (2008) in the case of Basel; Table 2], and a standard deviation of 10 per cent. Therefore, fluctuations to the stress field are introduced as observed from in situ stresses measurements (Valley & Evans 2009). In this way, each seed point is assigned a normally distributed random differential stress. This is one of the key stochastic elements of the model (Fig. 5a). Assuming a fracture orientation θ (i.e. the angle between σ 1 and the fracture normal) at each seed point (e.g. optimal or random orientation; Table 2), the seed points can be represented in a Mohr-Coulomb diagram with a normal and a shear stress component (Fig. 5b). Furthermore, a Mohr-Coulomb failure criterion expressed by cohesion c and a friction coefficient μ f is assumed (Table 2). During generation of the seed point population, stress states above the failure limit are rejected and redrawn. Different from Goertz-Allmann & Wiemer (2013), we here also define a criticality threshold that accounts for the fact that seed points cannot be too close to the failure limit. Stress states for which small changes on the order of 0.01 MPa are required for failure to occur, are likely to be triggered by small stress fluctuations, such as tidal forcing or passing seismic waves. We thus introduce a coefficient μ c , and replace (i.e. randomly redraw) seed points with stress states above a failure limit with friction μ f − μ c . This introduces a gap, here termed criticality, between the seed points closest to the failure limit and the failure limit itself (Fig. 5g). Criticality depends on normal stress and is about 0.55 MPa for normal stress of 55 MPa and μ c = 0.01. μ c is treated as free parameter to be calibrated against observations.
In a next step, a b value is assigned to each seed point using a negative linear relationship between b values and differential stress.
The assumption is based on observations in laboratory experiments (Scholz 1968;Amitrano 2003Amitrano , 2012, and from earthquakes magnitude distributions in different tectonic regimes (Schorlemmer et al. 2005;Gulia & Wiemer 2010) and at different depths (Gerstenberger et al. 2001;Spada et al. 2013). It reflects the observation that large magnitude events are more likely to occur at higher differential stresses. Goertz-Allmann & Wiemer (2013) also demonstrated that the assumption is necessary to reproduce observations at the Basel seismicity, such as the spatial distribution of b values as observed by Bachmann et al. (2012), and the high probability of larger events after shut-in. Once an event is triggered due to a pressure-induced change in effective stress, a magnitude is drawn from 10 5 random magnitudes forming a frequency magnitude distribution (FMD) with the corresponding b value at the triggered seed point. Randomly drawing magnitudes is another key stochastic element of the model. The relationship between b value and differential stress is defined by b max , the b value at zero differential stress. Note that such high b values are never reached as some differential stress is always needed for triggering an event. Following Goertz-Allmann & Wiemer (2013), the b value at large differential stresses (here σ > 135 MPa) is set to the ambient tectonic b value of 1.0 (Wiemer et al. 2009;Fig. 5h). Note that we use a magnitude of completeness M c = 0.9 in accordance with the conservative Figure 5. Explanation of stochastic seed model. (a) Random spatial distribution of seed points coloured according to the assigned differential stress. (b) b values at the same seed points using the relationship between b and σ (show in g). (c) For illustration a magnitude is assigned to all seed points. They are drawn from a catalogue of 10 000 events and a FMD with the corresponding b value. The second row Figs (d-f) shows the seed point in a, b and c in Mohr-Coulomb space. We assume optimal orientation for all seed points. (g) Inset figure of (d) showing that gap between the failure criterion and the seed point stresses created by the coefficient μ c . (h) Linear relationship between differential stress σ and b value. We adjust b max , although keeping the b value constant at 1.0 for differential stress higher than 135 MPa. value for the Basel seismic sequence defined by Bachmann et al. (2011), that is we only model events with magnitudes larger than 0.9. The value generally depends on the sensitivity of the monitoring network.
After failure at a seed point has occurred, a new stress state is assigned to the seed point using a drop in shear stress proportional to its normal stress; the proportionality constant is μ τ . For μ τ = 0.05 the stress drop would be τ = 2.75 MPa for a normal stress of 55 MPa, which is consistent with the average value of stress drop observed for events at Basel (Goertz-Allmann et al. 2011). Thus, multiple triggering of seed points is possible.

Parameter sensitivity
While a comprehensive study of the GSSM parameter sensitivities is beyond the scope of this paper, we shortly summarize the main effect of the b max , μ τ , μ c , μ f , c and N s : (1) The choice of b max does not affect the total number of events that occur in the model. However, it exerts first-order control on the overall b value of the modelled event sequence, and hence on the probability of larger events to occur. The lower b max the more likely large events occur.
(2) μ τ defines how easily events are re-triggered. Small stress drops arising from a small value of μ τ promotes re-triggering of seed points. As differential stress decreases progressively with each repeated triggering, the overall b value increases and thus the probability for larger events decreases slightly for smaller μ τ .
(3) The criticality coefficient μ c has an important effect on both the temporal evolution of seismicity and the overall b value. Because smaller overpressure is required to trigger events at low μ c , the postinjection pressure diffusion can trigger more events. Therefore, it has a major impact on seismicity after shut-in.
(4) μ f and c have similar effects: mean differential stress of all seed points is lower for lower cohesion and friction. Thus, the overall b value increases and fewer large events occur.
(5) Finally, N s affect the model results in a trivial manner: the higher N s the more events are triggered.
Note that there is a trade-off in estimating b max , μ f or c from seismicity observations as they all mainly have an effect on the overall b value of a sequence but not on the number of events. Estimating friction and cohesion from in situ observations is problematic, and a rough assumption has to be made if they are required as model input parameters. Later, during model calibration, we will adhere to an a priori assumption of friction and cohesion and only adjust b max to fit seismicity observations.

M O D E L C A L I B R AT I O N S T R AT E G Y
This section describes the procedure for calibrating our model against available observations of hydraulic behaviour and induced seismicity during stimulation experiments. The decoupled nature of the model makes it possible to calibrate both the flow and the seismicity models independently. In a first step, we use a pre-stimulation hydraulic test to estimate initial hydraulic properties. We calibrate permeability by minimizing the rms value of the difference between observed and modelled pre-test pressure values. We compute a number of possible parameter combinations of r eff and S, that fit the pre-stimulation pressure curve equally well. We thus keep all possible models and later derive the best-fitting combination using the extent of the seismicity cloud in the next step. In a second step, we use the stimulation pressure curve to fit the parameters C u , u t , p t and p t . We use a Monte Carlo simulation over the free parameter space to find the parameter sets that well reproduces the observations. In a third step, we use the calibrated pressure models in combination with the GSSM model to reproduce the following two characteristics of observed seismicity: (1) the temporal evolution of seismicity including the FMDs at all times, and (2) the spatial extent of the seismicity cloud. Regarding the latter, we do not attempt reproducing geometric details of the seismicity cloud, because this cannot be achieved with a radial-symmetric flow model.
We need to define a suitable metric to evaluate the performance of the modelled seismicity catalogues in comparison to the Basel observation. Following the experience gained within the Collaboratory of the Study of Earthquake Predictability framework (CSEP; Zechar et al. (2010), we use the log-likelihood L metric. We assume that the number of events with magnitudes within a certain magnitude bin and a certain time interval follows a Poisson distribution and estimate an expected rate λ ij from all model realizations (i and j are the number of the magnitude and time bins, respectively). For each bin, we then compute the likelihood L ij that the observed number of events can occur at a given modelled rate λ ij and calculate the total log-likelihood L as the sum over the logarithm of L ij over all bins. L is a negative number, which is larger the more the observed earthquake catalogues resembles the modelled realizations of seismicity catalogues. We need at least 500 model realizations to obtain a stable log-likelihood value.
With the log-likelihood L as measure for fit quality, we use again a Monte Carlo simulation to find a parameter set that reasonably reproduces the observations. We randomly vary b max , μ τ and μ c . We computed models for 300 parameter sets and for each parameter set 500 model realizations. The maximum extent of the modelled seismicity cloud is then compared to the observed one, that is, we compare the distance, which encloses 95 per cent of all events. In case of a misfit between modelled and observed extents, we adjust the combination of r eff and S accordingly, rerun the pressure model and redo the seismicity calibration procedure until the extent of the modelled seismicity cloud fits well the observed one. Thus, the seismicity cloud is fit only by adjusting S and r eff .

C A S E S T U DY: B A S E L
The model described in this paper is applied to the data set observed at the Basel EGS project in 2006 (Häring et al. 2008). The Basel EGS project aimed to become one of the first commercial power plants based on deep geothermal heat extraction from crystalline rock. It was planned to enhance reservoir permeability at about 4-5 km depth by injecting fluid at high pressure over a time period of more than 2 weeks. A seismic monitoring system was installed along with a hazard and risk management scheme-the 'traffic light system' suggested by Bommer et al. (2006). The monitoring system included six borehole sensors at depths of 300-2700 m depth. More than 900 events with magnitude larger than the magnitude of completeness (M c = 0.9) were recorded and located (Bachmann et al. 2011).
Injection rate was increased in a stepwise manner, until maximum injection rates of 57 l/s were reached on the fifth day of stimulation (Fig. 6c). Shortly after, an event of magnitude M L 2.6 occurred, which led to reduction of the injection rate, and, a few hours later, to the total stop of injection. A M L = 3.4 (corresponds to M w = 3.2) event, widely felt in the city of Basel, occurred about 5 hr later. The aversion of population and media against the project caused by these earthquakes led to temporal suspension of the experiment. In 2009, the project was fully cancelled as a consequence of a comprehensive risk study (SERIANEX, Baisch et al. 2009a). Allegedly, damage caused by the earthquakes included mostly fine cracks in plaster, which corresponds to an EMS intensity V. Insurance claims by homeowners reached about 7 million CHF most of which were also paid for.
In the following subsections, we demonstrate that our model approach can reasonably reproduce recordings of well-head pressure along with the injection rate, as well as the observed seismicity during the Basel stimulation experiment. We also show the model performance as a near-real-time forecasting tool with a pseudoprospective comparison between modelled and observed event rates. Note that we use the M w magnitude scale hereafter.

Injection rate and pressure
We obtain initial hydraulic model properties for Basel by using a pre-stimulation test performed several days before the actual stimulation, as also done by Ortiz et al. (2011). Fig. 6(a) shows the three-step pre-stimulation injection schedule. Fig. 6(b) compares the modelled and observed wellhead pressure during pre-stimulation injection. Prior to the start of the pre-stimulation test, the observed wellhead pressure showed an offset from zero pressure and a linear increase of 4.7 Pa s −1 (Baisch et al. 2009a). Thus, the observed pressure curve was set to zero at the beginning of injection and corrected for the linear increase. As suggested by Ortiz et al. (2011) for the Basel well, we use a wellbore storage coefficient of C ws = 9 × 10 −8 m 3 Pa −1 , which corresponds to an average compressibility of 4.1 × 10 −10 Pa −1 of 220 m 3 water residing in a 5000 m borehole. It accounts for the fact that the fluid column in the borehole can be compressed, which has a stronger effect the larger the fluid volume is (i.e. the longer the borehole is; Economides & Nolte 2000). The best-fitting initial permeability is κ 0 = 6.61 × 10 −18 m 2 ( Table 1). The value represents the equivalent porous medium permeability averaged over the open-hole section of h = 371 m. The corresponding average transmissibility is 2.45 × 10 −15 m 3 , which is on the same order of magnitude as the value ∼5 × 10 −15 m 3 reported by Häring et al. (2008). A series of best-fitting models were calculated for r eff between 1 and 5 m, while keeping κ 0 constant at 6.61 × 10 −18 m 2 . The combination that best fits the extent of the seismic cloud (about 600 m radius) is r eff = 2 m and S = 5.14 × 10 −12 Pa −1 (note that this is later derived from spatial extent seismicity). The difference between the modelled and the observed pressure curve is less than 0.3 MPa, with the largest deviation occurring during the first injection step. Note that Ortiz et al. (2011) suggest that the flow geometry may be best characterized by bilinear flow, based on the first step of the pre-stimulation test. Thus, our deviations during the first step may arise from a partly bilinear flow regime governing flow at early times, whereas radial flow, as assumed in our model, offers a suitable approximation after the first step.
Figs 6(c) and (d) show both the observed and the modelled stimulation pressure curve. Best-fitting parameters are given in Table 1. Note that u t = u t, = 135 , which implies that u can grow up to 270. In our model, a maximum value of about 230 is reached, (i.e. permeability has increased by that much). Häring et al. (2008) report a stimulation factor of 400 by applying standard well-test analysis of both pre-and post-stimulation tests. Although their factor is higher than the best-fitting value found here, they are in the same order of magnitude. The values for p t = 8 MPa and p t = 3.5 MPa found here imply that permeability enhancement emerges at 4.5 MPa, and reaches a maximum at 11.5 MPa. Fig. 6(f) shows the first 12 hr of stimulations and compares non-linear and linear flow models (i.e. with and without permeability enhancement). It shows that permeability enhancement becomes effective at about 7 MPa. The model fits the observations reasonably well given the strong assumption of a 2-D homogenous radial-symmetric continuum. Deviations between modelled and observed pressures are smaller than 4 MPa. The model fit is particularly weak after shut-in of the well, where the observed pressure decays gradually to pre-stimulation pressure, whereas the model pressure decreases much faster and stays at a nearly constant level of about 5 MPa. Fig. 7 shows the seismicity modelled with the GSS model, using the best-fitting parameters in Table 2, along with estimates of the  Fig. 7(d) shows the temporal evolution of the b value. We stress here that the increase of b values during injection and the decrease after shut-in could only be reproduced through the assumption of a relationship between the local seed point b value and differential stress. Fig. 7(e) shows the injection rate for comparison. Fig. 7(f) shows the observed and modelled FMD after 12 d. The modelled and observed b values of the sequence are comparable, although the observed one is slightly lower. Fig. 7(g) shows that the modelled and observed extents of the seismicity cloud are in agreement.

Seismicity
We demonstrate the model forecasting performance with a pseudo-prospective approach similar as Bachmann et al. (2011) did for statistical seismicity models. In a pseudo-prospective test of a seismicity model only a limited time interval of data is used to calibrate free model parameters. The calibrated model is then used to forecast the subsequent observed seismicity. By comparing the forecasted modelled data and the observed data, the forecasting capability of the model can be evaluated (e.g. Woessner et al. 2010). In our case, both the flow model and the GSS model were calibrated against only the first 48 and 72 hr of observations (referred to as 'learning periods'), and used to forecast the wellhead pressure (Fig. 8a) and the seismicity (Figs 8b and c) for the remaining time of the 12 d period. The modelled injection pressure for a learning period of 48 hr overestimates the maximum injection pressure, which is less severe for a learning period of 72 hr (Fig. 8a). Rates of seismicity (i.e. number of events larger than M c = 0.9), are also overestimated for the 48-hr learning period (Fig. 8b). For the 72-hr learning period the seismicity rates are well reproduced. The model performance does not improve substantially for longer learning periods.
In Fig. 9, we show the consistency between model and observations as a function of time and magnitude. Fig. 9(a) shows the injection rate as time reference. Fig. 9(b) indicates, if the observed cumulative number of events above different magnitude levels (like shown in Figs 7a-c) is consistent within the range given by all model realizations (dark grey; i.e. within the 95 per cent CI). The model is inconsistent on the first day, where it underestimates the number of events below M ≤ 1.6 (black colours). Similarly, the model underestimates magnitudes between 1.6 and 2.6 on the days 3-5 (white colours). Note that shut-in occurred at 5.7 d. The same colour scheme is used Fig. 9(c), where event numbers are grouped in magnitude bins of 0.1 and time bins of 12 hr. In most of the bins models and observations are consistent. However, patches of inconsistency occur; most of them represent underestimates by the model and only a few are overestimates. Fig. 9(d) demonstrates the consistency between model and observations similar as in Fig. 9(c) but for the pseudo-prospective test. The learning periods cover 1-5 d with steps of 1 d. For the learning periods of 1-4 d the subsequent 24 hr were forecasted and tested against observations. For the 5-d learning period, the entire remaining time period was forecasted. After the 2 d learning period, the forecasts are consistent with observations for most magnitude bins.

S E I S M I C H A Z A R D A N A LY S I S
Following Bachmann et al. (2011) and Mena et al. (2013), we can convert seismicity rates into time-dependent PSHA, expressed here as probability of exceeding a certain ground motion intensity. Rates in each time interval are converted using standard hazard calculation (Cornell 1968), using the EMS intensity attenuation relationship by Faeh et al. (2011). Results can be shown either as hazard curves for a given time interval, or as the probability of reaching or exceeding a given intensity.
One of the challenges when working with stochastic methods is to ensure that they well represent also rare and very rare events, because these may be the most relevant for hazard and risk assessment. One possibility is to extrapolate the simulated FMD from the numerous small events to the assumed maximum possible magnitude M max . However, as we show in the comparison between an extrapolation using the b value and the FMDs of 500, 1000 and 4000 model realizations (Fig. 10a), these two approaches deviate substantially for magnitudes above 3. The simulated FMD's are not adhering strictly to a power law, an effect already observed by Goertz-Allmann & Wiemer (2013) that is related to the difference in pre-and post-stimulation b values. The resulting hazard curves also vary substantially (Fig. 10b), with higher hazard in the intensity range IV-V because of the higher rates of magnitude 3-4 events. Beyond intensity VI and magnitude 4, the realizations are not stable enough because too few events are simulated. Fig. 10(c) shows a histogram of the maximum magnitude that occurred in each realization. The median of this distribution is M = 2.8, and can be regarded as the expectation value for the maximum observed magnitude during the 12 days during and after the Basel injection. The actual maximum observed magnitude was M w = 3.2.  Figure). Forecasting periods are 1 d except for the last learning period for which the entire remainder of the sequence was forecasted.
As can be seen from the cumulative distribution of maximum observed magnitudes in Fig. 10(d), it corresponds to a 20 per cent chance that a M w = 3.2 or larger event occurred in Basel. Based on a 95 per cent CI, a maximum observed magnitude between 2.4 and 3.9 would have been consistent with our model results.

A LT E R N AT I V E I N J E C T I O N S C E N A R I O S
Using the best-calibrated model, we now explore different injection scenarios (in terms of injection volume, injection pressure and depth of the reservoir) for their implication in terms of induced seismicity and resulting seismic hazard. Variation in total injected fluid volume is explored using the original Basel injection curve but terminating the injection as soon as 20-80 per cent (in steps of 20 per cent) of the maximum injected volume V max at Basel is reached (Fig. 11a). We also compute a scenario in which the total injected volume is pumped out of the well instead of only reducing the injection rate. Although such a scenario may not be applicable in reality, it attempts to explore the question if extracting fluid can help limiting seismic hazard. The results indicate that a slight reduction in hazard is achieved when pumping out water (Fig. 11e). Reducing the injected total fluid volume, however, has a substantial effect on the number of events (Figs 11b and c). The b value of the sequences does not vary, although the productivity (a value) changes proportional to the injected volume (Fig. 11d). This has a large impact on expected ground motion intensity (Fig. 11e).
We next compute scenarios with a constant total volume, but different injection pressure and rates. We arbitrarily choose a total volume corresponding to 60 per cent of V max at Basel; one scenario uses exactly the injection rate at Basel stopped at 60 per cent V max . In addition, we compute four scenarios with each four injection steps. We stretched and squeezed their injection duration so that it lasts 1/2, 3/4, 1.5 and 2 as long as the Basel 60 per cent V max scenario, and also adjusted the flow rate to obtain always the same total injection volume (Fig. 12a). The simulations show that varying injection pressure mainly affects the number of smaller events: high injection pressure allows for triggering more seed points that are not close to failure and tend to have lower differential stress. Thus, smaller events become more likely, whereas the number of larger events remains about the same (Figs 12b-d). Note that the total number of events in case of the highest pressure scenario is ∼60 per cent higher than in case of the lowest pressure scenario (Fig. 12d). However, the effect on the seismic hazard of the entire sequence is marginal (Fig. 12e). We conclude that different injection rates that all lead to the same total injected volume temporarily lead to different seismic hazard during injection, but eventually leads to the same cumulative seismic hazard.
We finally explore the effect of injection depth by varying σ 1 , σ 3 and hydrostatic pressure using the linear stress profile by Häring et al. (2008), to mimic stress states at 2500, 3000, 3500 and 4000 m depth. Average differential stress thus decreases from 110 to 60 MPa for 4500 and 2500 m, respectively. σ 3 changes from 75 to 42 MPa. The effect on the b value, and thus on the probability for large events and the hazard curves (Fig. 13e), is very strong, whereas the productivity only increases marginally (Figs 13b and d). Note that at 3000 and 2500 m depth, σ 3 is exceeded by the modelled pressure at some seed points, implying that tensile fracturing can occur. Tensile failure of pre-existing fractures is generally thought to be very inefficient in radiating seismic energy, and thus have only very small magnitudes if any (likely less than 0.9). Such events are removed from the synthetic catalogues. They contribute only 1 and 5 per cent to all events for 3000 and 2500 m, respectively. Thus, the number of events M ≥ 1 does not systematically increase for lower depth but decreases for depth of 3000 and 2500 m depth (Fig. 13b).

Effect of non-linear flow on modelling-induced seismicity
Sensitivity analysis of both the flow and the seismicity model components indicate that the spatial extent of the seismicity cloud induced by fluid injection depends not only on the initial hydraulic properties. The pressure disturbance penetrates farther into the reservoir for lower stimulation pressures p t (Figs 4e and f). Thus, the final size of the seismicity cloud depends on the minimum pressure required for stimulation to become effective. Similarly, the criticality in the stochastic seed model strongly affects the extent of the seismicity cloud. The closer the stress state at a seed point lies to the failure limit, the smaller the pressure required for triggering. Thus, decreasing μ c has two effects: (1) more events can be triggered after shut-in, and (2) the spatial extent of seismicity increases.
In our case, we fit the spatial extent of seismicity only by adjusting the specific storage coefficient S, because this is the most sensitive parameter regarding spatial extent. The best-fitting value of 5.14 × 10 −12 Pa −1 is low for any rock mass. It requires a rock mass bulk modulus of up to 200 GPa, if it is derived by applying linear storage using the bulk moduli of rock and fluid (e.g. Rutqvist & Stephansson 2003). A possible reason for such a low value of S may be that viscosity η is set too high (see eq. 1). In the relevant pressure regime of 45-75 MPa, a value of η = 2.5 × 10 −4 Pa s is realistic for temperature of 115-120 • C, whereas the value may decrease to η = 1.5 × 10 −4 Pa s for temperatures as high as 190 • C (Likhachev 2003) corresponding to the undisturbed temperature of the reservoir (Häring et al. 2008). We thus expect that viscosity may be somewhat lower than 2.5 × 10 −4 Pa s. However, this effect cannot fully explain the low value of S. Other parameters relevant for defining the spatial extent of seismicity, such as p t and μ c , may play a role in altering the estimate of S.
Nevertheless, we can use the estimated value of S to calculate diffusivity of the reservoir along with viscosity η and permeability κ and obtain a value for diffusivity of D 0 = 0.0051 m 2 s −1 before, and D stim = 1.18 m 2 s −1 after stimulation. Both values strongly deviate from the effective diffusivity D eff = 0.05 m 2 s −1 estimated for the Basel reservoir by Dinske (2010). Even for higher, that is more realistic, values of S, the obtained unstimulated and stimulated diffusivity values are not in agreement with D eff . As the D eff lies in between those two diffusivity values, it is not clear what it actually represents. We conclude that estimating hydraulic parameters from observed seismicity cannot be based on linear spherically symmetric diffusion as suggested by Shapiro et al. (1997Shapiro et al. ( , 2002. For modelling flow in a stimulation context the full set of physical processes relevant in induced seismicity need to be considered comprehensively (Bruel 2007;McClure & Horne 2011); at a minimum, however, one needs to consider non-linear irreversible pressure diffusion.

Implication of model simplifications
By formulating non-linear pressure diffusion with a pressure dependent permeability, we enhance computational efficiency as required for the use in PSHA and real-time applications, but neglect thermohydromechanical processes that constitute the underlying physics of induced seismicity. In reality, permeability is enhanced through slipinduced shear dilation that is triggered by fluid pressure (Rutqvist & Stephansson 2003). Although our strategy of using irreversibly pressure-dependent permeability is able to reproduce wellhead pressure reasonably well, it remains uncertain whether the actual spatial distribution of permeability enhancement is realistic as well. As a consequence, interpretation of the model scenarios  is limited with respect to the reservoir permeability. We may draw conclusions on how to limit seismicity using different injection strategies based on those scenarios; however, the models may not predict if a certain reduction of seismic hazard can still produce the volume and permeability enhancement needed in order to achieve the required reservoir properties. Our current model therefore does not fulfil the third requirement for future reservoir simulators outlined in the introduction (Section 1.3).
Another limitation of the model is introduced by not explicitly considering stress redistribution after slip at a seed point, which affects the stress state at other seed points. As shown by Catalli et al. (2013) and Schoenball et al. (2012), stress redistribution after induced earthquakes, usually termed Coulomb failure stress change ( CFS; Stein 1999), is relevant for understanding induced seismicity. Values of CFS are in the order 0.1-1 MPa in the vicinity (at <100 m distance) of hypocentres in the case of reservoir stimulation. Although the stress magnitude of such changes is lower than the fluid pressure that induces events, it plays a role in modifying the stress state of nearby potential earthquake sources. Catalli et al. (2013) showed that CFS becomes more important as stimulation progresses and induced slip accumulates, and is most relevant for post-injection events. Our model produces slightly too low event rates around and after shut-in (Fig. 9c), which may well be explicable by the shortcoming of ignoring stress interaction between seed points. Further development of our modelling approach may have to account for such effects of stress interaction, although the use in real-time application would require hypocen-tre precision of better than 50 m (Catalli et al. 2013), which are challenging.

Scaling of seismicity with injection parameters
The results in Fig. 12 indicate that seismic hazard depends only marginally on the injection pressure. Different injection histories at constant total volume but varying rate mainly results in higher numbers of smaller events (M ≥ 1), whereas larger events (M ≥ 3) are about equally probable to occur (Fig. 12). For high-pressure injections (i.e. high injection rates), overpressure is especially high close to the injection point, whereas at larger distances the pressure is only marginally higher compared to low-pressure injections. Therefore, more events are triggered close to the injection point, and they tend to have smaller magnitudes as they have smaller differential stress. However, the maximum extent of the seismic cloud does not change as long as volume is kept constant. Although for different injection rates and pressures the seismic hazard varies temporarily, the total seismic hazard (i.e. the hazard covering the injection period and the pressure relaxation time) is independent on the injection pressure and rate. The results are in agreement with the ones by McClure & Horne (2011), which suggest that higher injection pressure results in higher seismicity rates but not in larger magnitudes. Their scenarios did not keep total injection volume constant while varying injection pressure. Thus, in their model, the higher seismicity resulting from higher injection pressure also includes a volume effect that we also observe in our model. Our models also suggest that post-shut-in events cannot be limited by slowly reducing injection pressure instead of abruptly shutting-in. In contrast, it may be advantageous to stop injection immediately and even allow venting of the well as was done in the case of Basel, because that limits the fluid volume that enters the reservoir-a strategy also suggested by McClure & Horne (2011).
In our model, the seismic hazard strongly depends on the total injected volume (Fig. 11). We further illustrate this dependence of seismic hazard on injected volume in Fig. 14. The black line shows the mean and the gray shadings the 95 and 99 per cent CI of the maximum magnitude of each model realization (as in Fig. 10c). To cover a wider volume range in Fig. 11, we added more scenarios, that is we computed two Basel-type scenarios where we kept injecting at maximum injection rate instead of shutting-in, as well as four scenarios in which only 10, 5, 3 and 1 per cent of V max at Basel was injected. The expected maximum observed magnitude scales linearly with the logarithm of injected volume. The results support the volume dependence of seismic hazard that is also assumed in the Shapiro model (Shapiro et al. 2010), and helps explaining why it performs well in forecasting seismicity (Mena et al. 2013).
We also show in Fig. 14 the maximum observed magnitudes of injection experiments summarized by Evans et al. (2012) for European cases, and added prominent beyond-Europe cases: Cooper Basin, Australia (Baisch et al. 2006(Baisch et al. , 2009b, Paradox Valley, USA (Ake et al. 2005) and Ogachi, Japan (Shapiro et al. 2007). All cases are represented with different colours to show different injection depths, and with different marker shapes to indicate the tectonic seismic hazard at each site (i.e. peak ground acceleration PGA that is expected to be exceeded in 50 yr with a probability of 10 per cent). Considerable scatter of the various case studies conceals any possible trend of these maximum observed magnitudes with volume. Variable injection depth and tectonic setting may explain some of the scatter. Most case studies that lie around the trend lines defined by the Basel-calibrated model have an injection depth of 4 km or greater. Exceptions are the cases Cesano and Torre Alfina, which are shallower, but are located in an area with moderate tectonic seismic hazard as in Basel. In contrast, all cases that lie far below the Basel-calibrated model trend lines, have either low tectonic seismic hazard or were injected at shallow depth. In the outstanding case of Ogachi, injection was performed at 1 km and at a site with high tectonic seismic hazard. The depth dependence of the induced seismicity is an important input for the design of future experiments. The strong dependence of seismicity on depth that our model scenarios suggest (Fig. 11) is an immediate consequence of the implemented relationship between b values and differential stress. Fig. 15 summarizes the effect on the maximum expected magnitude, and on the seismic hazard for depths varying between 1 and 5 km. It illustrates that not only the expected maximum magnitude increases strongly with depth, but also the uncertainty thereof (Fig. 15a). Interestingly, also the expected ground motion intensity at the surface increases with injection depth, illustrating that the increase with depth in the number of larger events outweighs the attenuation effect. However, we stress that the strong dependence of simulated seismicity with depth in our model is a function of the assumed relationship between differential stress and b values, which is only calibrated for the Basel injection. It is possibly not well constrained outside of this depth range; the calibration with additional case studies from a range of depth before drawing definitive conclusions are required. Nevertheless, we argue that a depth dependence of seismicity is reasonable, because (1) a b value on depth dependence is seen in tectonic seismicity Figure 14. Scaling between injection volume and the expected maximum observed magnitude. Also shown is the 95 and 99 per cent CI. The maximum expected magnitude scales linearly with the logarithm of injected volume. Also shown are case studies from Evans et al. (2012), and following additional ones: Ogachi, Japan (Shapiro et al. 2007), Cooper Basin, Australia (Baisch et al. 2006(Baisch et al. , 2009a; Paradox Valley, USA (Ake et al. 2005). The markers are coloured according to injection depth and have shapes corresponding to the ambient seismic hazard expressed as the PGA with 10 per cent probability of being exceeded in 50 yr. Further, we show the modelled maximum expected magnitude for depths of 3.5 and 2.5 km (Fig. 15a) for comparison.  Fig. 10c). Also shown is the 95 per cent and the 99 per cent CI. (b) Depth dependence of hazard expressed as the EMS intensity exceeded with a probability of 99, 10 and 1 per cent. The shorter travel distance of seismic waves at shallower depths is outweighed by the lower magnitudes expected due to a depth dependence of b values. (Gerstenberger et al. 2001;Spada et al. 2013) and (2) also the a value, the average tectonic activity, decreases towards shallower depths (Spada et al. 2013).
We acknowledge that also various other site specific parameters determine the level of induced seismicity. For tectonic events a dependence of b values on the tectonic stress regime was reported (Schorlemmer et al. 2005;Gulia & Wiemer 2010). A measure for tectonic predisposition (like stress regime or stressing rates), and other site specific factors was suggested by Shapiro et al. (2010) by introducing the seismogenic index. It expresses the readiness of a reservoir to react by seismic energy release. In our models, the seed density largely defines the number of events that can occur, whereas the relationship between seed point b values and differential stress determines the occurrence rate of large events. Both these model parameters represent site-specific predisposition similar to the seismogenic index. Our model results and also the strong scattering of observed maximum magnitudes shown in Fig. 14 confirm that seismicity not only scales with volume (McGarr 1976), but strongly depends on the potential energy from tectonic stressing stored in the system before stimulation.
Although site specific conditions may explain some of the scatter of observed magnitudes in Fig. 14, we suggest that already the inherent stochastic nature of events induced in a critically stressed crust explain much of the scatter. This scatter is best seen in the histogram of Fig. 10(c), where we show the distribution of expected largest event magnitudes for 4000 Basel simulations. In these simulations, all injection parameters are kept constant (depth, volume, pressure, b values, etc.) and the only reason for the variability between individual simulations is due to the randomly drawn distribution of faults and their sizes. The Basel experiment could have resulted with equal probability in a maximum magnitude of 2.6 rather than the observed 3.2, or any magnitude between 2.4 and 3.9 that are within the 95 per cent CI. This wide range of possible maximum observed magnitudes arises purely from the stochastic nature of the distribution of potential failure planes in a critically stressed Earth crust.
We believe that it is important to clearly distinguish the 'maximum observed magnitude' M(obs) max , for one experiment at a certain site, and the 'maximum possible earthquake', M max , as typically used in PSHA studies (Wiemer et al. 2009). M(obs) max represents the mean of the distribution of many experiments or simulations; in other words, what is, on average, the maximum event that is expected for each simulation. In the case of Basel, the value is 2.8 (Fig. 10c). The second one represents the very end of the upper tail of this distribution, which is an extremely rare event that may happen only every 1000, or 10 000 times. In the case of our simulation, we assume as M max the regionally assumed M max of 7.2 (Wiemer et al. 2009), but 4000 simulation only once gave a value as large as 4.5. To sample on average a M max of 7.2 once, we would need to draw more than 1 million times, and such rare events trend to be not relevant for estimating the hazard at commonly used probability levels. As it is the case of classical PSHA, M max is a parameter that is very difficult to estimate based on observed data (Coppersmith et al. 2009), because observation periods are generally much too short to sample M max . It is therefore possible that physical mechanisms not included in our model would truncate the magnitude distribution below the typically assumed tectonic M max values. However, we argue that the observed maximum magnitudes for the about 20 cases of injection into crystalline shown in Fig. 14 carry no information on M max . Thus, using compilations of case studies as in Fig. 14 to define the upper truncation of the FMD and to be used in seismic hazard studies (see also McGarr (1976); SERIANEX by Baisch et al. 2009a) is incorrect and may lead to a substantial underestimation of the seismic hazard. Nevertheless, such plots, especially when combined with modelling results, are important because they provide information on the expected M(obs) max for a given injection volume.
The same conclusions can be drawn when arguing from a physicsbased rather than a statistical point of view: In a critically stressed crust, the maximum possible moment release during an injection cannot be computed deterministically from the moment release required to accommodate volumetric expansion of the reservoir, as suggested first by McGarr (1976). Instead, a finite probability remains that induced slip occurs at a highly critically stressed nucleation patch that may dynamically grow into an event that has a larger moment release than predicted by the injected volume alone. The statement is also in agreement with the study by McClure & Horne (2011) and by Garagash & Germanovic (2012). The latter show conceptually that quasi-static slip along a pressurized fault plane can grow into a dynamic slip in areas that are not limited by the pressurized area. Their results indicate that the maximum possible magnitude is not defined by the pressurized volume alone, but instead also depends on initial stress conditions. As those parameters always remain largely unknown, induced seismicity hazard should be discussed within a probabilistic instead of a deterministic framework.

S U M M A RY
The most essential conclusions of our study are summarized in the following: (1) We present a non-linear flow model that accounts for an irreversible pressure-dependent increase of permeability. The model is able to reproduce the observed wellhead pressure at Basel using the employed injection rate. The modelled reservoir pressure responsible for inducing seismic events is at a realistic order of magnitude, which is not possible with linear pressure diffusion models.
(2) The generation of numerous synthetic earthquake catalogues including event magnitudes is accomplished with a stochastic process. At seed points that are triggered by overpressure, we randomly draw magnitudes from FMDs with b values that depend on a predefined differential stress at the seed point. We thus make use of the observation that b values of tectonic earthquake catalogues reveal a dependency on differential stress.
(3) The proposed two-step modelling approach is able to reproduce the observed temporal and spatial evolution of seismicity rates at different magnitudes. The model can also reproduce the observation that seismic hazard remains high after shut-in to some degree. Also temporally and spatially variable b values, as observed by Bachmann et al. (2011Bachmann et al. ( , 2012) arise from our model.
(4) Alternative stimulation scenarios are produced using the model calibrated against Basel data. We find that seismicity scales strongly with injected volume and less so with injection pressure. The strongest effect arises from different injection depth. The result, however, strongly depends on the assumed relationship between b value and differential stress. Regarding design of injection strategies of future EGS, we recommend limiting both injection depth and volume. At a given injection volume, short injections at high pressure may be an advantage as more small events relative to large events are induced. We further recommend stopping injection by immediately shutting-in and venting the well to reduce the fluid volume in the reservoir.
(5) We argue that scaling relationships cannot be based on single magnitudes (e.g. maximum observed magnitude), but rather on the probability of certain magnitudes to be exceeded. Deterministic estimates of the maximum observed magnitude-for example as a function of injected volume-are currently not applicable, because the complexity and unpredictability of tectonic predisposition and earthquake occurrence forces us to rely on probabilistic methods.
(6) Although physical considerations in our model are limited, it represents a first attempt to combine numerical reservoir modelling with PSHA -two fields of research that often diverge in their methods and philosophy. We stress, however, that coping with seismic hazard in reservoir engineering projects cannot work without merging those two fields.

A C K N O W L E D G E M E N T S
We are grateful to GeoPower Basel for permission to publish these results. We thank Andres Alcolea, Flaminia Catalli and Arnaud Mignan for useful comments on the manuscript. Also we thanks to Nicolas Hummel and an anonymous reviewer for their useful comments and suggestions. This work was supported by the European Commission through the funds provided within the FP7 project GEISER, grant agreement no. 241321-2.