Static behaviour of induced seismicity

The standard paradigm to describe seismicity induced by fluid injection is to apply nonlinear diffusion dynamics in a poroelastic medium. I show that the spatiotemporal behaviour and rate evolution of induced seismicity can, instead, be expressed by geometric operations on a static stress field produced by volume change at depth. I obtain laws similar in form to the ones derived from poroelasticity while requiring a lower description length. Although fluid flow is known to occur in the ground, it is not pertinent to the behaviour of induced seismicity. The proposed model is equivalent to the static stress model for tectonic foreshocks generated by the Non- Critical Precursory Accelerating Seismicity Theory. This study hence verifies the explanatory power of this theory outside of its original scope.


Introduction
Induced seismicity is a growing concern for the energy-sector industry relying on fluid injection in the deep parts of the Earth's crust [Ellsworth, 2013;Mignan et al., 2015]. At the same time, fluid injection sites provide natural laboratories to study the impact of increased fluid pressure on earthquake generation [Majer et al., 2007]. with n a positive integer [Shapiro and Dinske, 2009]. These two descriptive laws can be derived from the differential equations of poroelasticity [Biot, 1941] under various assumptions. The full description of the process requires complex numeric modelling coupling fluid flow, heat transport and geomechanics [Rutqvist, 2011]. These models, numerically cumbersome, can become intractable because of the sheer number of parameters [Miller, 2015]. Attempts to additionally correct for the known discrepancies between Biot's theory and rock experiments have led to a large variety of model assumptions [Berryman and Wang, 2001], indicating that poroelasticity results are ambiguous.
I will demonstrate that a simple static stress model can explain the two empirical laws of induced seismicity without requiring any concept of poroelasticity.
The proposed theoretical framework hence avoids the aforementioned shortcomings by suggesting an origin of induced seismicity that does not involve fluid flow in a porous medium (although fluid flow indeed occurs). Historically, a similar static stress model was proposed for the tectonic regime under the Non-Critical Precursory Accelerating Seismicity Theory (N-C PAST) [Mignan et al., 2007;Mignan, 2012]. Its application to induced seismicity data will allow a more fundamental investigation of the relationship between static stress and earthquake generation. To test the model, I will use data from the 2006 Basel Enhanced Geothermal System (EGS) stimulation experiment including the flow rate of injected fluids [Häring et al., 2008] and the relocated catalogue of induced seismicity [Kraft and Deichmann, 2014].

The Non-Critical Precursory Accelerating Seismicity Theory (N-C PAST)
The N-C PAST has been proposed to explain the precursory seismicity patterns observed before large earthquakes from geometric operations in the spatiotemporal stress field generated by constant tectonic stress accumulation [Mignan et al., 2007;Mignan, 2012]. In particular, it provides a mathematical expression of temporal power-laws without requiring local interactions between the elements of the system [Sammis and Sornette, 2002;Mignan, 2011]. Therefore earthquakes are considered passive (static) tracers of the stress accumulation process, in contrast with active earthquake cascading in a critical process (hence the term "non-critical"). The concept of self-organized criticality [Bak and Tang, 1989] is seldom used to explain induced seismicity [Grasso and Sornette, 1998]. Since there is no equivalent of a mainshock in induced seismicity, the criticality versus noncriticality debate has limited meaning in that case. However, the underlying process of static stress changes considered in the N-C PAST can be tested against the observed spatiotemporal behaviour of induced seismicity.
The N-C PAST postulates that earthquake activity can be categorized in three regimes -background, quiescence and activation -depending on the spatiotemporal stress field σ(r,t) defined from the boundary conditions σ(r " +∞, t) = σ 0 * and σ(r = 0, t) = σ 0 + t + σ 0 * , with h the depth of the fault segment base, r the distance along the stress field gradient from the fault's surface projection, σ 0 < 0 the stress drop associated to a hypothetical silent slip occurring at t 0 at the base of the fault, the tectonic stress rate on the fault, σ 0 * the crustal background stress, n = 3 the spatial diffusion exponent for static stress and t f the mainshock occurrence time [Mignan et al., 2007] (Fig. 1a).
Background, quiescence and activation regimes are defined by event densities δ b0 , δ bm , and δ bp for |σ| ≤ σ 0 * ± Δσ * , σ < σ 0 * -Δσ * and σ > σ 0 * + Δσ * , respectively, with the boundary layer ±Δσ * the background stress amplitude range. By definition, δ bm < δ b0 4 < δ bp with each seismicity regime assumed isotropic and homogeneous in space (i.e. role of fault network neglected). Correlation between earthquake productivity and static stress changes is well established [King, 2007]. The distinction of three unique seismicity regimes with constant event density, the main assumption of the N-C PAST, is discussed later on.
Figure 1: Seismicity spatiotemporal behaviour described by the N-C PAST static stress model (tectonic case [Mignan, 2012]): (a) Spatiotemporal evolution of the stress field σ(r,t) generated by constant stress accumulation on a fault located at r = 0 (Eq. 1). Background, quiescence and activation seismicity regimes are described by densities of events δ b0 , δ bm , and δ bp for |σ| ≤ σ 0 * ± Δσ * , σ < σ 0 * -Δσ * and σ > σ 0 * + Δσ * , respectively; (b) Temporal evolution of quiescence and activation envelopes r * (t) with σ(r * ) = σ 0 * ± Δσ * (Eq. 2); (c) Rate of events µ(t) in a disc of constant radius max(r*) (Eq. 3); (d) Cumulative number of events N(t) (Eq. 4) of power-law form In the tectonic case, static stress changes are underloading due to hypothetical precursory silent slip on the fault at t 0 followed by overloading due to hypothetical asperities delaying rupture on the fault after t p * [Mignan, 2012]. The three seismicity regimes are then defined as solid spatiotemporal objects with envelopes The parameters t m * = t mid -Δσ * / and t p * = t mid + Δσ * / represent the times of quiescence disappearance and of activation appearance, respectively, with σ(0, t mid ) = σ 0 * . The background seismicity regime is defined by subtracting the quiescence and activation envelopes r A * (t) and r Q * (t) from a larger constant envelope r max ≥ max(r * ) (Fig. 1b). While trivial along , concepts of geometric modelling may be required to represent these seismicity solids in threedimensional Euclidian space [Gallier, 1999] in which the vector is possibly curved [Mignan, 2011]. The non-stationary background seismicity rate µ(t) is then defined in the volume of maximum extent r max by with k a geometric parameter and d the spatial dimension. For the tectonic case in which r max >> h, the volume is assumed a cylinder with k = π, d = 2 and δ the density of epicentres in space (Fig. 1c). Finally, the cumulative number of events N(t) is defined as which represents a power-law time-to-failure equation of the form the first term representing the linear background seismicity and the second term the quiescence or activation power-law behaviour observed prior to some large mainshocks ( Fig. 1d) [Sammis and Sornette, 2002].

Application of the N-C PAST static stress model to induced seismicity
In the case of an EGS stimulation, the stress source is the fluid injected at depth with overpressure where K is the bulk modulus, ΔV the volume change per time unit and V 0 the infinitesimal volume subjected to pressure effect per time unit at the borehole located at r = 0. The injected volume V(t) is determined from the flow rate profile Q(t), as with t 0 the starting time of the injection. The change of volume is then defined as with Δt a time increment.
In the EGS case, r ≅ h with h the borehole depth and induced seismicity defined as hypocentres. The spatiotemporal stress field σ(r,t) becomes with r the distance along the stress field gradient from the borehole, n = 3 the spatial diffusion exponent for static stress and r 0 → 0 the infinitesimal radius of volume V 0 = kr 0 d /t 0 , t 0 = 1 being the time unit. Activation represents the case when fluids are injected and quiescence when fluids are ejected (bleed-off). It follows that which suggests that the spatiotemporal shape of the induced seismicity envelope depends on the nth-root of the flow rate profile Q(t) (with n = 3 in the static stress case). This parabolic relationship is similar to the generalized form r(t) ∝ m(t) 1/d derived from nonlinear poroelasticity in a heterogeneous medium where m is the cumulative mass of injected fluid and d the spatial dimension [Shapiro and Dinske, 2009]. The main difference between the two physical approaches is in the underlying stress field, which is here static and in poroelasticity, dynamic and related to the displacement gradient of the fluid mass [Rudnicki, 1986]. It is trivial to derive Eq.
The induced seismicity rate µ(t) is then defined by Eq.
(10), k = 4π/3 and d = 3, assuming a spherical spatial volume (i.e. isotropic stress field). For the activation phase (i.e. stimulation period), it follows that 8 The induced seismicity case d = n = 3 confirms the linear relationship between cumulative injected volume and cumulative number of induced earthquakes N(t) ∝ V(t) previously derived from poroelasticity [e.g., Shapiro and Dinske, 2009]. In contrast with poroelasticity, this second law is a direct consequence of the first. The d = n condition also yields the simplified form of Eq. (10) where the one free parameter is the normalized background stress amplitude range  . 2a). The N-C PAST thus predicts an activation envelope r A * for t 0 ≤ t < t 1 and a quiescence envelope r Q * for t ≥ t 1 (Eq. 13). The activation and quiescence envelopes are fitted to the Basel data using Δ * ∈ [10 -3 , 10 -1 ] day -1 (light curves) and Δt = 1/4 day. The results are shown in Figure 2b. The value Δ * = 0.007 day -1 (dark curves) provides the best fit to the data, defined from the best score S = (w A +w Q )/2 with w A and w Q the ratio of events of distance r ≤ r A * and r ≥ r Q * in the injection and bleedingoff phases, respectively. Figure 2c shows S as a function of Δ * for Δt = {1/12, 1/8, 1/4} day, which indicates that the results remain stable for lower time increments.  [Häring et al., 2008]; (b) Spatiotemporal distribution of relocated induced seismicity [Kraft and Deichmann, 2014] with r the distance from the borehole. The activation and quiescence envelopes r A * (t) and r Q * (t)

Application to the 2006 Basel EGS induced seismicity sequence
are defined from Eq. (13) with parameters Δ * = 0.007 day -1 (dark curves) and Δt = 1/4 day. The light curves represent the range Δ * ∈ [10 -3 , 10 -1 ] day -1 in 0.1 increments in the log 10 scale. Points represent the induced earthquakes, which colour indicates how they are declared; (c) Score S = (w A +w Q )/2 with w A and w Q the ratio of events of distance r ≤ r A * and r ≥ r Q * in the injection and bleeding-off phases, respectively. The vertical line represents Δ * = 0.007 day -1 .
I evaluate δ b0 = 10 -10 event/m 3 /day by counting all earthquakes declared in the national Swiss catalogue (ECOS-09 1 ) and located within 10 km of the borehole of coordinates (7.594°E; 47.586°N) and depth 4.36 km. It means that ~1 tectonic earthquake is expected in average in the space-time window considered. Due to the low tectonic activity in the area, I approximate δ b0 = δ bm = 0 event/m 3 /day (i.e., total quiescence). The theory shows a good agreement with the observations with 97% of the seismicity below r A * during the injection phase (red points in Fig. 2b) and 98% of the seismicity above r Q * during the bleeding-off phase (orange to yellow points).
The density of events above r Q * is however not δ b0 but which represents the temporal diffusion of induced seismicity with τ the average time constant. Eq. (14) represents a relaxation process from the overloading state to the background state. The results here suggest that only the events declared as background (grey points) and quiescence events (blue points) are outliers. The observed variations in r below r A * and above r Q * are not explained by the model, which only predicts the behaviour of the activation and quiescence fronts. The second-order variations may be due to anisotropic effects and for t > t(max(r A * )) to additional spatial diffusion effects. Figure 3 shows the 6-hour rate of induced seismicity µ(t) and the cumulative number of induced events N(t), observed and predicted. With δ b0 = δ bm = 0 and taking into account induced seismicity temporal diffusion, the rate of induced seismicity where δ bp = 4.68 10 -7 event/m 3 /day (production parameter) and τ = 1.18 day (diffusion parameter) are obtained by maximum-likelihood estimation (MLE), set S t = {Δt, …, Eq. (15) infers that induced seismicity is fully explained by overloading, in agreement with the observation of no causal relationships between events in the Basel sequence [Langenbruch et al., 2011]. The predicted rate (Eq. 15) and predicted cumulative number of events (Eq. 4) fit the data well, as shown in Figures 3a and 3b, respectively. The role of temporal diffusion is observed after t 1 -Δt and is the only contributor to induced seismicity after t 1 . Of three functional forms tested to describe diffusion (exponential, stretched exponential and power law), the exponential (Eq. 14) was verified to be the best model for the Basel case (following the formalism and tests proposed by Clauset et al. [2009]).

Conclusions
I have demonstrated that the two principal induced seismicity descriptive laws can be explained from simple geometric operations in a static stress field without requiring any concept derived from poroelasticity. The two descriptive laws had been previously obtained by considering the differential equations of poroelasticity [Biot, 12 1941;Rudnicki, 1986] under different assumptions [Shapiro and Dinske, 2009], which indicates that the static stress model defined from algebraic expressions requires a lower description length [Kolmogorov, 1965]. This is crudely inferred here from the difference between the lengths of the present demonstration and of published poroelasticity demonstrations [e.g., Shapiro and Dinske, 2009]. Histogram of the observed 6-hour induced seismicity rate µ(t) with fit based on Eq.
I also showed that the controlling parameter is the normalized background stress amplitude range Δ * , which questions the usefulness of permeability and diffusivity parameters in induced seismicity analyses and might explain why these parameters remain elusive [Miller, 2015]. In that view, permeability could depend on the "external loading configuration" instead of on the material itself, as recently proposed in the case of the static friction coefficient [Ben-David and Fineberg, 2013].
Testing of the model on other induced seismicity sequences will determine if Δ * is itself universal, region-specific or related to the static stress memory of the crust, hence if Δ * depends or not on the tectonic loading configuration at EGS natural laboratory sites. Similar questions apply to the earthquake production parameter δ bp and if the two parameters are independent or correlated.
The main assumption of the N-C PAST is to consider three unique seismicity regimes (quiescence, background and activation) defined by the event productions δ bm < δ b0 < δ bp . There are two possible physical alternatives to justify this choice: (1) it represents the fundamental behaviour of the Earth crust, which would hence act as a capacitor, with strain energy storage and δ bp analogues to electrical energy storage and capacitance, respectively; (2) the proposed step function is a simplification of the true stress-production profile, which remains unknown and is so far best characterized by three regimes [e.g., King, 2007]. Both alternatives allow defining spatiotemporal solids over which geometrical operations yield algebraic expressions of the induced seismicity behaviour.