Stress‐Based and Convolutional Forecasting of Injection‐Induced Seismicity: Application to the Otaniemi Geothermal Reservoir Stimulation

Induced seismicity observed during Enhanced Geothermal Stimulation at Otaniemi, Finland is modeled using both statistical and physical approaches. The physical model produces simulations closest to the observations when assuming rate‐and‐state friction for shear failure with diffusivity matching the pressure build‐up at the well‐head at onset of injections. Rate‐and‐state friction implies a time‐dependent earthquake nucleation process which is found to be essential in reproducing the spatial pattern of seismicity. This implies that permeability inferred from the expansion of the seismicity triggering front (Shapiro et al., 1997, https://doi.org/10.1111/j.1365-246x.1997.tb01215.x) can be biased. We suggest a heuristic method to account for this bias that is independent of the earthquake magnitude detection threshold. Our modeling suggests that the Omori law decay during injection shut‐ins results mainly from stress relaxation by pore pressure diffusion. During successive stimulations, seismicity should only be induced where the previous maximum of Coulomb stress changes is exceeded. This effect, commonly referred to as the Kaiser effect, is not clearly visible in the data from Otaniemi. The different injection locations at the various stimulation stages may have resulted in sufficiently different effective stress distributions that the effect was muted. We describe a statistical model whereby seismicity rate is estimated from convolution of the injection history with a kernel which approximates earthquake triggering by fluid diffusion. The statistical method has superior computational efficiency to the physical model and fits the observations as well as the physical model. This approach is applicable provided the Kaiser effect is not strong, as was the case in Otaniemi.


Injection-Induced Seismicity: Mechanisms and Forecasting Methods
Induced seismicity can result from either a stress or strength change on a fracture or fault. The effect of injection is generally assessed by considering pore pressure diffusion in the medium and the consequent decrease in the effective normal stress as according to Terzaghi's principle (Skempton, 1984). This first-order description of the stress state has been effective in explaining various aspects of induced seismicity, including the √ evolution of the seismicity front (Shapiro et al., 2006(Shapiro et al., , 1997 and general spatiotemporal patterns of induced seismicity (Elmar & Shapiro, 2002;Shapiro et al., 1999Shapiro et al., , 2002 as early as the pioneering study at the Rangely oil field (Raleigh et al., 1976). An additional step in the description of stress changes due to a fluid injection is the theory of poroelasticity which describes the coupling between fluid flow and deformation of the solid skeleton. Poroelasticity has been shown to play a role in triggering earthquakes in addition to pore pressure evolution (Segall, 1989;Segall et al., 1994;Segall & Lu, 2015), particularly outside the characteristic pore pressure diffusion length (Goebel & Brodsky, 2018;Zbinden et al., 2020). Although the magnitude of stress changes from poroelasticity is estimated to account for typically only about a tenth of that from pore pressure diffusion (Zhai & Shirazei, 2018), its consideration is often required for complete explanations of the observed seismicity in space and time.
A fluid injection can result in "hydrofractures" (Mode-I opening fractures) or shear fractures (Mode-II or Mode-III). Induced earthquakes generally result from shear failure. While linear elastic fracture mechanics is commonly employed in modeling the growth of cracks in Mode-I and the consequent stress changes, modeling shear failure requires an appropriate friction law. One kind of models is based on the Mohr-Coulomb failure criterion in which slip occurs once the ratio of the shear stress to the normal stress on a fault reaches a pre-defined threshold, the static friction coefficient, and drops to the dynamic friction coefficient either at the immediate onset of slip or gradually with fault slip. However, there is ample evidence from laboratory studies and natural observations that the initiation of slip involves in fact a gradual decrease of friction associated with aseismic slip, often referred to as the nucleation process. Such an evolution of friction is commonly described using the rate-and-state friction law derived from frictional sliding experiments in the laboratory (Ampuero & Rubin, 2008;Dieterich, 1994;Dieterich & Linker, 1992;Marone, 1998;Ruina, 1983).
The non-instantaneous nucleation process implied by rate-and-state friction can explain a number of phenomenological observations such as the Omori decay of seismicity rate during aftershocks (Dieterich, 1994) or the low sensitivity of seismicity to solid-earth tides (e.g., Beeler & Lockner, 2003). The rate-and-state formalism has also shown success in explaining the relationship between stress and seismicity rate due to diking (e.g., Toda et al., 2002) and aseismic slip (e.g., Segall et al., 2006). In the context of induced seismicity, rate-and-state friction has been applied to explain certain nonlinear features such as the time lag between induced seismicity and stress perturbations (e.g., Candela et al., 2019;Dempsey & Riffault, 2019;Norbeck & Rubinstein, 2018;Richter et al., 2020). It is important to note that, in principle, the activation of a fault by a pore pressure increase does not necessarily imply seismic slip (e.g., Guglielmi et al., 2015). In fact, there is observational evidence that injection-induced fault slip is mostly conditionally stable (Bourouis & Bernard, 2007;Calò et al., 2011;Goodfellow et al., 2015;Guglielmi et al., 2015;Scotti & Cornet, 1994), as is expected from the nucleation model based on rate-and-state friction and that seismicity is in fact occurring outside the zones of high pore pressure (Cappa et al., 2019;De Barros et al., 2018;Wei et al., 2015).
More specifically with regards to hydraulic stimulation of geothermal wells, important questions arise regarding the differences between the Mohr-Coulomb and rate-and-state friction-based models considering the rapid stressing rate that is common in such operations. Mohr-Coulomb models coupled with linear slip weakening can result in realistic simulations of seismic ruptures while accounting for the nucleation process (Olsen et al., 1997). This is not the case for single-degree-of-freedom spring-slider systems often employed for modeling induced seismicity. The commonly used model of Dieterich (1994) based on rate-and-state friction can converge to models based on the Mohr-Coulomb criterion at the rapid equilibrium limit. It is also possible that rate-and-state effects on nucleation may be significant at the relatively short timescale of intense injection cycles during stimulation.
A hysteresis effect, often referred to as the Kaiser effect, is also commonly observed in induced seismicity. The Kaiser effect refers to the observation when a material submitted to a series of loading cycles of increasing amplitude fails gradually, further failure generally occurs at a stress level exceeding the maximum stress reached in previous cycles. This effect explains the observation that acoustic emissions during rock failure stop if the stress decreases and do not resume until the medium is loaded to its previous maximum (Lavrov, 2003). How a nucleation source "remembers" its loading history has proven to be essential in reproducing various observations in induced seismicity, such as time delays of the seismicity rate in response to perturbations of the injection rate and regions of seismic quiescence behind triggering fronts (Baisch et al., 2010(Baisch et al., , 2006Dempsey & Riffault, 2019).
Numerous physical models have been developed to incorporate stress changes, pore-pressure changes, and failure mechanisms in a single framework (Gaucher et al., 2015;Grigoli et al., 2017). A notable example of physical models that accounts for rate-and-state friction in particular, is presented by Segall and Lu (2015), where changes in stresses by fluid injections into an infinite poro-elastic medium were used as input to the model of Dieterich (1994), relating seismicity and stress rates among a population of nucleation sources. Although the framework was originally used to investigate poroelastic effects during shut-in and to address the common observation that maximum magnitude events often occur after injections cease (Grigoli et al., 2018;Häring et al., 2008), it can be used more generally to study induced seismicity in response to various injection scenarios (e.g., Zhai & Shirazei, 2018). Finite-fault and fracture network models accounting for rate-and-state friction have also been developed (Almakari et al., 2019;Dublanchet, 2018;Larochelle et al., 2021;McClure & Horne, 2011) to examine rupture properties and the effect of heterogeneous fault properties on the seismicity rate. Numerous factors make it difficult, however, to resort to such models in practice, such as the high computational cost of solvers and poor resolution of pre-existing heterogeneities in the sub-surface-in particular, the distribution of stress and strength-with a level of detail that cannot be constrained with observation. Some representations of heterogeneities are essential in reproducing well-established statistical properties of earthquakes (Dempsey et al., 2016;Zoller et al., 2005) such as the Gutenberg-Richter law which describes the magnitude-frequency distribution of earthquakes (Gutenberg & Richter, 1956).
Due to the complexity of stress-based models along with the difficulty to calibrate the model parameters, a number of studies have alternatively explored data-driven statistical modeling. Such models often hinge on the Gutenberg-Richter law (Gutenberg & Richter, 1956) and the assumption that earthquakes follow a Poisson process. Additionally, they often model earthquake triggering as a cascading process based on the Omori law (Utsu, 2002) which fits commonly observed patterns of the decay of seismicity rate during aftershock sequences. A popular example is the epidemic type aftershock model (ETAS; e.g., Ogata, 1988), which represents the total seismicity as a linear superposition of homogeneous Poisson processes, to represent mainshock and aftershock sequences (e.g., Bachmann et al., 2011;Lei et al., 2008;Mena et al., 2013). Such models have the advantage of resulting in very realistic synthetic catalogs since they incorporate statistical properties directly derived from observations. However, statistical approaches are in principle less transportable from one reservoir to another as they lack explicit connections to the mechanical and hydro-geological properties of the medium. The development of hybrid models that account for the complex network of physical mechanisms while being generalizable and applicable to various injection sites and scenarios is therefore an active area of research (Gaucher et al., 2015).

Data Presentation and Analysis
The seismic catalog analyzed in this study comes from a geothermal well stimulation project operated by St1 Deep Heat Ltd. near the campus of Aalto University in Otaniemi, Finland and is compiled by Leonhardt et al. (2021). The injection well (OTN-3 in Figure 1) was drilled to a depth of 6.1 km into Precambrian crystalline (gneiss and granite) rocks. Approximately 18,000 m 3 of water was injected over the course of 49 days from 4 June to 22 July 4 of 27 in 2018. The injection history was divided into five successive stages moving upward from the bottom of the well (Figure 1). Pumping parameters of the injection such as the injection rate and well-head pressure were tuned as part of a Traffic Light System (TLS), the details of which are presented in Ader et al. (2020) and Kwiatek et al. (2019). The stimulation consisted of numerous cycles of injections and pauses of varying duration. The injection history also included periods of bleed-off's where injection was stopped and backflow out of the well was allowed.
The stimulations were monitored with surface and borehole seismometers providing excellent detection and location of the induced earthquakes (Hillers et al., 2020;Kwiatek et al., 2019). Namely, the monitoring network consisted of a seismometer array at 2.20-2.65 km depth in a separate well (OTN-2), located around 400 m from OTN-3, in addition to a 12-station network installed in 0.3-1.15 km deep wells (Figure 1). The catalog consists of 61,150 events in total ( Figure 2) and 1986 relocated events with spatial uncertainty of ±52 m (Figure 3). The magnitude of completeness is estimated to be M c = −1.1.
A few salient features of the observed seismicity guide our modeling. First, the seismicity rate has a positive correlation to the injection rate in time, accompanied by finite periods over which it increases and decreases in response to injections and shut-ins, respectively. We indeed note that the seismicity rate reaches a similar magnitude for injections far apart in time but equal in the flow rate. Second, the decay pattern in the seismicity rate, R, during injection pauses is well-matched by the Omori law where t is time, t r is the time it takes for the seismicity rate to halve, and R 0 is the seismicity rate at the onset of decay. A fit to one of the injection pause periods is shown in Figure 4. Note that the more general "modified Omori law" (Utsu, 2002) allows a 1/t p decay of seismicity rate; here the p-value is close to 1. The close match to the Omori law is consistent with observations of the decay rate in induced seismicity following shutins reported in a number of previous studies (Almakari et al., 2019;Bachmann et al., 2012Bachmann et al., , 2011Langenbruch & Shapiro, 2010). Lastly, the relocated catalog ( Figure 3) shows a rather diffuse distribution of seismicity,   , with c tf = c horner = 0.018 m 2 /s. It is difficult to assess a level of agreement between the triggering front and the relocated catalog given the limited sample size. Yet, clusters of events far beyond the curve suggest poroelastic triggering. It is also possible that they are due to leaks in the casing, as evidenced by their locations close to the well path shown in the vertical section view (bottom-left). In the map (bottom-right) and and vertical section views, the well is drawn in black with stimulated sections of the well and occurrence time of events color-coded correspondingly. M HEL refers to the local Helsinki magnitude scale. The color-coding reveals little correlation in space between events and stimulation stages.
suggesting that the injection stimulated fractures were distributed within a relatively large volume (∼ 1 km 3 ) around the open sections of the well by diffusion of pore pressure.
The exact origin of Omori law decay remains poorly understood; it could be due to the finite nucleation process governed by rate-and-state friction (Dieterich, 1994) or by instantaneous nucleation and postseismic creep that predict a p-value of approximately 1 (Perfettini & Avouac, 2004). This process was suggested to have occurred during a 10 MPa stimulation of a geothermal well at ∼ 3 km depth at Soultz-sous-Forêt (Bourouis & Bernard, 2007). Similarly, stress relaxation by pore pressure diffusion (Nur & Booker, 1972) predicts a seismicity decay also closely resembling the Omori law with a p-value typically between 1 and 2 (Langenbruch & Shapiro, 2010;Miller, 2020). Studying the properties of the Omori-like decay provides a valuable opportunity to re-examine its mechanical origins and the physical mechanisms that drive induced seismicity.

Linear Transfer Function and Convolution Model
The direct relationship between the injection and observed seismicity rate suggests that it may be represented by a linear transfer function of the injection history (Avouac et al., 2020). To quantify this relationship, we use the algorithm of Marsan and Lengline (2008) which was originally designed to determine the kernels characterizing how earthquakes trigger other earthquakes. The algorithm estimates weights as a function of distance and time which, after normalization, represent the probability that any earthquake was triggered by any previous earthquake. We adapted the algorithm here to determine the weight relating earthquakes to injections as the source of trigger. As justified later on, secondary triggering is ignored (i.e., aftershocks of triggered events are ignored). We assume that the observed seismicity rate density, λ(x, t), or the number of earthquakes in unit time can be modeled by a linear superposition of the influence from all previous injections such that: where λ 0 is the uniform background rate density, and λ i (t) represents the rate density at time t incurred by injection i. A nonlinear behavior may in reality arise from the possible coupling between fluid pressure and permeability, and from the seismicity model. Rate-and-state friction and the Kaiser effect are indeed sources of nonlinearity, as we discuss in greater detail below.
The kernel λ(Δt) (referred to as the bare rates) that defines λ i (t) is found through an iterative process: First, we begin with an initial guess for λ(Δt) and compute the triggering weights between injection i and event j, w i,j = α j λ(t j − t i ) and the background weight w 0,j = α j λ 0 where α j is a normalization coefficient to satisfy that ∑ −1 =0 = 1 . Here, w i,j = 0 if t i > t j (earthquakes cannot be triggered by future injections). Second, λ(Δt) is updated as follows: where A is the set of pairs such that |t j − t i | ≤ δt, and N is the number of total earthquakes. Thus, δt becomes the discretization parameter of the algorithm. The two main assumptions of the model are linearity of the rate density that allows superposition of λ i and the existence of a mean-field response to injections that is independent of event magnitude or injection volume. Demonstration of the algorithm on a simple synthetic catalog and its sensitivity to discretization parameters are illustrated in Text S1 in Supporting Information S1.
Injections are divided into individual cycles by binning them into regular 10 min intervals. The result reveals a time decay proportional to 1/t ( Figure 5). This is consistent with the observed Omori law decay following shutins and also with the period of build-up in seismicity at the beginning of injections. It is also possible to use this approach to estimate spatial kernels. The results are not presented here as we found the size of the data set and the quality of the locations to be insufficient to get well constrained kernels.  Figure 2). A short-period prior to shut-in is shown with a sky blue background. The shut-in period is indicated with a gray background. The decay pattern in seismicity rate during the shut-in is matched well with an Omori decay function (modified Omor-Utsu law with p = 1), plotted in light purple. The dotted lines and shaded areas in-between indicate the 95% confidence interval of the fit. The fitted value of t r and the bounds of the confidence interval of the fit are indicated in the legend.

of 27
The observation that the response to step-like decrease of injection rate leads to a 1/t Omori law decay can be used to estimate a Green's function, g(t) (Avouac et al., 2020). Since the derivative of a step function is a Dirac delta function, g(t) can be found by simply differentiating the Omori law in time The predicted seismicity rate can then be obtained from a simple convolution where R and u are the seismicity and injection rate, respectively. Bleed-off's are implemented as negative injection rates (likewise to all forthcoming models in this study). To construct the kernel for the specific case of Otaniemi, t r is chosen by fitting the Omori law to the last of the injection pauses of durations significantly longer than the average injection cycle (about 20 hr). Then, R 0 is determined so as to yield a total number of events equal to the number of earthquakes in the catalog. t r and R 0 are found to be 24.1 hr and 208.9 events per hour, respectively. Although Avouac et al. (2020) reported that the data suggests a systematic increase of t r during the stimulation likely due to the increasing volume of the domain of increased pore pressure, we use a constant value of t r as the resulting difference to the fit is minor.
The model result is displayed with the observed catalog in Figure 6a. It follows remarkably well the observed seismicity rate variations; bulk of the observed seismicity is included within the 95% confidence interval, calculated by assuming events are governed by a non-homogeneous Poisson process following the modeled seismicity rate. The model also closely matches the decay rate during injection pauses and the build-up rate at the onset of injection cycles.
To quantify the goodness of fit, we use both the Kolmogorov-Smirnov test (Massey, 1951) and the Poisson log-likelihood (Dempsey & Suckale, 2017). The Kolmogorov-Smirnov test returns the KS-statistic, which is the maximum difference between the cumulative distribution functions given by the prediction and the observation. The Poisson log-likelihood is the appropriate metric if earthquakes are assumed to result from a Poisson process, even if inhomogeneous in the case the rate varies in time and space. So the metric is valid as long as secondary aftershocks can be ignored. This assumption is tested by analyzing the distribution of interevent distances in space and time using the method of Zaliapin and Ben-Zion (2013). The result is shown in Figure S4 in Supporting Information S1, which displays a uni-modal distribution instead of the bi-modal distribution that would be expected in case of clustering due to aftershock sequences. This is consistent with the analysis by Kwiatek et al. (2019) which shows that aftershocks account for no more than 10% of the events in their seismicity catalog and the observation that aftershock sequences are rarely observed in seismicity induced by hydraulic stimulations (e.g., Baisch & Harjes, 2003). One advantage of the Poisson log-likelihood and the Kolmogorov-Smirnov test is also that the metrics do not require binning of the point process (Dempsey & Suckale, 2017). Binning is used in the figures only for convenience to represent the data. The log-likelihood function is given by where θ is the set of model parameters and t j is the occurrence time of event j = {1, 2, …, n}. We report the KS-statistic here, preferred to the log-likelihood which is sensitive to the choice of units for R, but we see good qualitative agreement between the two measures as summarized in Table 2. The KS-statistic for the convolution model returns 0.036. The quality of the fit is impressive considering the simplicity of the model-which involves only two parameters. It also contradicts the premise that various nonlinear mechanisms driving induced seismicity, such as the nonlinear ity of rate-and-state friction, the Kaiser effect, and changes in permeability due to high pore pressure and the development of hydraulic fractures, should result in a nonlinear response overall. It 9 of 27 may be that nonlinear effects in Otaniemi are in fact small despite the relatively large stress variations induced by hydraulic stimulation, the possibility of which we explore with our physical models later on and in Supporting Information S1.

Physical Modeling
We now present a physical model based on stress evolution from pore pressure diffusion and poroelasticity along with shear failure criterion following rate-and-state friction. The medium is treated to be infinite, homogeneous, and isotropic. Neglecting the effect of the free surface is justified by the relatively large depth of the injections compared to the dimensions of the seismicity cloud ( Figure 3). The induced stresses can then be calculated using the analytical solutions for a point source from Rudnicki (1986).
where p and σ ij are the pore pressure and stress tensor, and r and t the distance from injection source and time, respectively; λ u = 2μν u /(1 − −2ν u ) is the undrained Lamé parameter and the drained Lamé parameter without the subscript u; c is the hydraulic diffusivity which depends on permeability, k and viscosity, η. Here, we add the subscript "true" to k and c to distinguish between the true and apparent values of the diffusivity, the notions of which are explored in greater detail by our following analysis. We assume the point source is a good approximation of the injections in Otaniemi given the length of the stimulated wells relative to the size of the total stimulated volume. The model is nearly identical to that introduced by Segall and Lu (2015). Poroelastic properties which lack constraints from the field, along with a fixed fault-orientation are chosen as those in Segall and Lu (2015) to represent a general case. Ambient normal stress of 155 MPa is approximated using the average depth of the injection. All fixed parameters and their dimensions are listed in Table 1.
Stress changes become the input to the ODE formulation of Dieterich (1994), to solve for seismicity rate in space and time. The alternative integral formulation of Heimisson and Segall (2018) is used here as it is more tractable numerically for injection scenarios such as in Otaniemi that consist of abrupt onsets and shut-ins of injections While the global fit to the observations are comparable to other models, it lacks rapid variations of the seismicity rate in-between injection cycles compared to the rate-and-state models-evident of qualitative differences in modeling the stress state relative to failure and delayed nucleation mechanisms. All models (besides (c)) consistently capture temporal trends of the seismicity rate, such as the Omori-law decay during shut-ins and build-up periods at the onset of injections, with the linear convolution model requiring the fewest parameters and lowest computational cost. Model parameters and goodness-of-fit metrics are summarized in Table 2.
where r b is the background seismicity rate, the background stressing rate, a the rate-and-state friction parameter, σ the normal stress, 0 and τ 0 the initial effective normal and shear stress, and and τ the applied effective normal and shear stress, respectively. Synthetic catalogs are produced by sampling events from a non-homogeneous Poisson process using the acceptance-rejection method.
The Kaiser effect is inherent in the formulation of Dieterich (1994) and Heimisson and Segall (2018). This results from the fact that the nucleation process is delayed if the stress decreases and resumes once the stress gets back to its previous peak level. The Kaiser effect is clearly demonstrated if we use the model to compute the response of the seismicity rate to a sinusoidal stressing history ( Figure S5 in Supporting Information S1). The different injection locations must stimulate new volumes of rock and lead to new hydraulic pathways. So we might expect the Kaiser effect to be significant within a single stage but to be less relevant from one stage to the other. The impact of the Kaiser effect may be more appropriately represented by resetting the stressing history at the onset of each stage. To this effect, we start a new simulation with the same initial conditions and compound the results for the final catalog. This model is hereafter referred to as the rate-and-state model. Note that the validity of resetting the stress history could be questioned given that the seismicity clouds during the different stages largely overlap ( Figure 3) suggesting overlapping stimulated volumes.
We use the measured flow rates and pressure to estimate hydraulic diffusivity. An estimate of the diffusivity that fits the rate of pressure decay during injection pauses is made by the Horner analysis. Since the analytical solutions of the present model are derived for spherical flow in a 3-D medium, the conventional Horner analysis originally derived for 2-D flow into a vertically confined aquifer (Horne, 1995;Zimmermann, 2018) is adapted to be consistent with Equations 7 and 8. Details on the adaptation and fitting process are presented in Text S2 in Supporting Information S1. The analysis gives a diffusivity of c horner = 0.018 m 2 /s, and a global fit to the entire pressure history using a Gaussian likelihood function gives an effective well radius and ambient pore pressure of 44 m and 43.5 MPa, respectively. The model fits the measured pressure history well during the entire stimulation, especially during the injection pauses ( Figure 7a). A fit to the pressure history with diffusivity as a free parameter, however, gives a higher value of c bu = 0.044 m 2 /s (subscript "bu" standing for "build-up") that better matches the rate of pressure build-up at the onset of injection cycles (Figure 7b) with an effective radius and ambient pore pressure of 31 m and 54.9 MPa, respectively. c bu also over predicts the rate of pressure decay during injection pauses. While constraints on the effective radius-a measure of the damage zone surrounding the well that causes pressure drops-are difficult to quantify, ambient pore pressure in either cases are close to its bounds considering the temperature-dependence of fluid density at injection depth. When comparing the theoretical triggering front derived by Shapiro et al. (1997), that is, = √ 4 where c tf is the diffusivity chosen to draw the triggering front, c horner appears to fit the spatial extent of near-field events better (Figure 3). We therefore use c horner = c true as a starting point for the models and refer to its theoretical triggering front as the "reference triggering front." We revise this assumption later and note that the diffusivity derived from the Horner analysis fits the pressure drop at shut-ins, as should be the case by design, but does not match the pressure build-up when injections start again (Figure 7a).
The posterior distribution on the set of parameters associated to the seismicity model , and r b is found using the affine invariant Markov chain Monte Carlo (MCMC) Ensemble sampler of Goodman and Weare (2010) maximizing the log-likelihood given by Equation 6. In order to simplify the sampling process, the sampler computes the posterior of a and given that r b -which is a simple multiplicative factor to the normalized seismicity rate-is adjusted for each pair of a and to match the total number of observed events (61,150 events). The sampler conducts 2,000-5,000 iterations of 32 walkers with the chain length made to be longer 50 times the auto-correlation length in order to ensure full exploration of the posterior distribution. The prior is assumed to be uniform for both variables between the range of 10 −5 -10 −2 and 0.1-5 kPa/yr for a and , respectively, although the shape of the prior is seen to have little effect on the posterior given the large sample size.
, and r b of maximum likelihood is found to be 0.0002, 3.05 kPa/yr and 12.1 events/days, respectively, and the resulting model is shown in Figure 6b. The model follows the observations quite well in time, with a KS-statistic of 0.029, slightly lower than the value of 0.036 obtained with the convolution model. The model succeeds in reproducing the main temporal features of the observed catalog: (a) direct correlation between the injection and seismicity rate and (b) Omori-law decay during shut-ins. In space, the fit is much less compelling (Figure 8b). The triggering front lags significantly behind the reference triggering front with a much smaller mean of the distribution. Yet in both time and space, resetting of the stress history at each injection stage turns out to be essential in reproducing important features of the observation. The best fit using the model without resetting of the stress history (a = 0.0001, = 4.89 kPa/yr, and r b = 25.9 events/day) as shown in Figure 6c has relatively minimal seismicity rate during the second half of the injection history due to the Kaiser effect. In space, it is completely devoid of any seismicity close to the injection well during this period (Figure 8c). Far-field seismicity much beyond the reference triggering front is largely attributed to background stressing as poroelastic stress perturbations are small relative to pore pressure changes.

Adjusting Model Diffusivity to Spatio-Temporal Distribution of Seismicity
Given that the rate-and-state model fails to match the observations in space assuming the diffusivity inferred from Horner analysis, we now examine the possible underestimation of the diffusivity by the Horner analysis. Following the seminal study of Shapiro et al. (1997), it has become common practice to infer the diffusivity from fitting = √ 4 to the propagation of the seismicity front, or the triggering front-defined by the outline of the outermost events of the seismicity cloud extending from the well. However, we note that c tf of the rate-and-state model shows a significant mismatch by a factor of ∼3 from c true = c horner prescribed in the model (Figure 8b). This discrepancy is due to the role of delayed nucleation represented by aσ. As shown by Wenzel (2017), the parameter aσ of the rate-and-state model acts as a threshold triggering stress that restricts the extent of the triggering front. The sensitivity of the triggering front to aσ is clearly visible in Figure 9 which compares two synthetic catalogs that only differ in the prescribed values of a. In the scope of the rate-and-state model or stress thresholds as commonly used in Mohr-Coulomb models, inference of the diffusivity from the apparent migration of seismicity requires considerations of both c and a. Additionally, the method of inferring the diffusivity from the triggering front may depend on the earthquake detection thresholds. A higher detection threshold may give a more poorly resolved catalog in space that could lead to a different estimation of the triggering front. Furthermore, the position of the triggering front can be obscured even more by background seismicity and far-field events triggered by poroelastic effects. Fitting the seismicity front represented by the envelope of the seismicity cloud, places a lot of weight on potentially biased and not particularly well-defined features.
In consideration of such complications, one would wish for a definition of the seismicity front that is independent of the number of events in the catalog and robust to factors of discrepancy between observations and model predictions. We therefore propose an approach to infer c true from the spatial distribution of the seismicity as opposed to the triggering front. A simple way is to fit the distribution as a function of distance and time from the point of injection with a known analytical expression. We recall that the half-norm distribution is the solution to the diffusion equation in response to a Dirac point source in a 3-D medium where the standard deviation of the distribution, Λ(t), is a function of time such that: This inspires our approach to fit Equation 10 to the rate-and-state model in response to a constant injection scenario. The half-norm distribution indeed turns out to provide a relatively good fit ( Figure 10); it matches well the bulk of the distribution but tends to slightly overestimate seismicity rate at larger distances. Indeed, we do not make the claim that the half-norm distribution is the best possible fit and acknowledge there may be other distributions that could better match the rate-and-state model although they are not explored further here. Furthermore, plotting the evolution of Λ versus time reveals that it follows closely √ . We make the assumption that the remaining discrepancy can be modeled as a multiplicative factor such that: where {l} is a set of non-dimensional parameters. Thus, c hg is a measure of the radial spreading of the seismicity relative to the point of injection ("hg" standing for half-Gaussian distribution). In order to apply this method to Otaniemi, we attempt to estimate c hg from the relocated catalog. One disadvantage of the method is that it requires a set of relocated events large enough to constrain the evolution of c hg with confidence. As detailed in Text S3 in Supporting Information S1, we can indirectly estimate from the cumulative relocated catalog giving c hg = 0.011 m 2 /s ( Figure S6 in Supporting Information S1).
We find the relationship γ h (l) empirically by observing the systematic dependence of γ h on l as reproduced by the rate-and-state model. We assume l depends not only on pore fluid transport properties but also rate-and-state properties such as a. We find to be relevant the ratio l = aσ/p q , where = 4 0 is the characteristic pore pressure for given injection rate q, and L is the size of the computational domain. Higher values of aσ would produce a stronger threshold effect and suppress seismicity migration, the extent of which would depend on its strength relative to the induced pressure, p q . A series of single boxcar injections are simulated for a range of c and a. We find a rational function of aσ/p q that fits γ h as shown in Figure 11. Although the reason for the exact functional form of the relationship is not obvious, the quality of the fit is compelling. The observed trend is also consistent with the known role of aσ: Higher values of a suppresses seismicity at further distances, decreasing c hg and consequentially γ h . The functional fit allows new uncertainty estimates of the diffusivity in Otaniemi. Figure 11 shows the difference between the predicted and true values of diffusivity for a range of c true and a, given the estimated value of c hg = 0.011 m 2 /s and an injection rate, q = 10 L/min. Although this is a much lower injection rate than the average in Otaniemi there are also significant differences between the idealized boxcar injections used to produce Figure 11 and the much more complex schedule in Otaniemi. One can see that accounting for the role of delayed nucleation significantly widens the possible range of diffusivity in Otaniemi. Namely, the functional fit considers The fit to space in the rate-and-state model is significantly improved with c true = c bu = 0.044 m 2 /s. The rate-and-state models consist of far-field seismic activity, although mostly from background stressing distributed uniformly in space rather than through a systematic variation from poroelastic stress perturbations. (e) The Coulomb model with c true = c bu = 0.044 m 2 /s significantly overpredicts the distribution of seismicity in space as does the theoretical triggering front for c tf = c bu , suggesting that the role of delayed nucleation on seismicity migration is essential in reproducing the observed spatio-temporal evolution of seismicity in Otaniemi given the likely diffusivities. Model parameters and goodness-of-fit metrics are summarized in Table 2. Figure 9. Sensitivity of triggering front to delayed nucleation: Synthetic catalogs for two parameter sets only differing by a (0.0001 and 0.001 in top and bottom, respectively) are shown. Lower a, which translates to lower aσ, results in a much further extent of the triggering front, due to the role of delayed nucleation that acts proportionally to a threshold stress for the triggering of events as explained in detail by (Wenzel, 2017). Along with the reference triggering front in red, an additional √ 4 curve is drawn in orange for a = 0.001, with c tf modified by a factor of 0.3 that better matches the apparent triggering front. equally likely much higher values of c true than would be predicted by the triggering front observed in Otaniemi given sufficient rate-and-state effects.
In light of this finding, we test the possibility that c bu = 0.044 m 2 /s is in fact closer to c true in Otaniemi than c horner as the inconsistency between the triggering front using c bu = c tf and the relocated catalog are borne due to rate-and-state effects. We test this hypothesis by finding the best fit of the rate-and-state model using c bu = c true . The effective radius and ambient pore pressure are adjusted to 31.1 m and 54.9 MPa, respectively, to best fit the well pressure measurements. The resulting fit for the seismicity rate in time is shown in Figure 6d, and the corresponding synthetic catalog in space is shown in Figure 8d. a, , and r b are found to be 0.00006, 1.29 kPa/yr, and 4.7 events/day, respectively. The fit in time bears no significant improvement from the fit using c horner = c true , although the KS-statistic is slightly lower at 0.025. The fit in space is much improved with a higher mean of the distribution and cluster of events that encompasses greater portions of the relocated catalog. One region the model performs rather poorly on is capturing the back-propagation front starting around the 500 hr mark. Itis possible that the back-propagation fronts, whose occurrence in time would correspond to the drawdown periods used for the Horner analysis, is still governed by the lower diffusivity c horner . It could be that the back-propagation consists of two separate migration patterns, based on the observation that the initial portions of the back-propagation front Figure 10. Evolution of spatial distribution of seismicity for rate-and-state model: Spatial profiles of the seismicity rate are plotted in blue at various times for the rate-and-state model in response to a single boxcar injection. Half-norm distributions, in green, are used to fit the model-generated distribution. The line style is alternated between solid and dashed between each time step for clarity. The half-norm distributions evolve with a time-dependent shape parameter, Λ(t), which closely follows √ as shown in the inset of the top figure. Figure 11. Inference of diffusivity accounting for role of delayed nucleation on seismicity migration: An empirical relationship for the multiplicative factor, γ h , of Λ(t) = √ ℎ is found in terms of the non-dimensional ratio aσ/p q (left). The fit can be used to infer new uncertainty estimates on the diffusivity of the medium given apparent spreading of the radial distribution of the seismicity in Otaniemi, that is, c hg = 0.011 m 2 /s. Contour plot on the right shows the percent difference between the true diffusivity and the predicted diffusivity from the functional fit γ h (aσ/p q ) for a range of a and c true . Considerations of the role of delayed nucleation on seismicity migration makes higher diffusivities more likely than previously considering solely the theoretical triggering front of Shapiro et al. (1997). are predicted quite well by the model (starting at around the 450 hr mark). This could be due to a propagation of the seismicity governed by different mechanisms than pore pressure diffusion, such as stress transfer by aseismic slip (Dublanchet & De Barros, 2021), although it is difficult to constrain the exact mechanism of seismicity migration given their possibly similar characteristics ( ∼ √ ) .
The differences between c bu and c horner may be indications of distinct hydraulic processes that govern the well-head pressure and the spatial distribution of seismicity. One could imagine that the well-head pressure is more representative of the diffusivity of the medium immediately surrounding the well. On the other hand, the spatial distribution of seismicity may be more dependent on the path of highest hydraulic conductivity within the entire stimulated volume. The abrupt cessation of seismic activity close to the injection well following shut-in could be associated to a decrease in the diffusivity due to fracture healing, leading to the lower estimate of c horner . It is also important to note that the two diffusivities require different values of a, , and r b , such that their independent measurements would provide stricter constraints on c true . We see that the higher estimate c bu inferred from this analysis yields synthetic catalogs in better agreement with the observed seismicity in time and space. We conclude using the triggering front to infer the diffusivity may yield a significantly biased estimate if the effect of earthquake nucleation is ignored.

Design of the Spatio-Temporal Convolution Kernel
We now use the physical model as a basis to extend the temporal convolution model to space. We look for a new kernel with spatial dependence such that the convolution is as follows: The spatial component of the kernel is constructed by using the half-norm distribution, as identified in Section 6, with a shape parameter increasing as √ ℎ . Combining with the Omori law as the temporal component as previously gives the integral of the kernel which is differentiated in time to obtain the response to Dirac forcing The three parameters of the model are c hg = 0.011 m 2 /s, R 0 = 213.5 events/hr, and t r = 28.5 hr, as estimated from the data. The fit to the temporal evolution of seismicity is, by design, identical to the fit obtained with the kernel in time introduced earlier (Figure 6a). The model provides now in addition a good match to the observations in space, especially with regards to the triggering and back-propagation fronts (Figure 8a). Overall, the convolution method approximates the physical model and fit the observations quite well, albeit with a drastically shorter computing time-by at least an order of magnitude-thanks to the use of the fast Fourier transform (the convolution is transformed into a simple product in the Fourier domain).

Comparisons of Coulomb and Rate-And-State Models
Both rate-and-state and Mohr-Coulomb models are widely used in modeling induced seismicity. The standard Coulomb model assumes a population of faults with a uniform distribution of initial stress up to the maximum shear stress allowed by static friction (e.g., Ader et al., 2014). We show in supplement that this simplest version of the Coulomb model does not fit the observations neither in time nor in space (Text S4 and Figure S7 in Supporting Information S1). A number of studies which have tested the applicability of the Coulomb model to induced seismicity found it necessary to introduce a stress threshold that needs to be exceeded for earthquake triggering (e.g., Bourne et al., 2018;Dempsey & Riffault, 2019;Dempsey & Suckale, 2017;Langenbruch & Shapiro, 2010;Rothert & Shapiro, 2003). The physical justification for the inclusion of the threshold, hereafter referred to as C cpt , is to account for the population of faults activated during the stimulation that were initially in a relaxed state of stress, not close to failure. In this case, triggering would be delayed due to their initial strength excess rather than due to the nucleation process. The explanation is probably relevant in stable tectonic areas (e.g., Bourne et al., 2018;Dempsey & Riffault, 2019;Dempsey & Suckale, 2017;Langenbruch & Shapiro, 2010). Wenzel (2017) demonstrates the response of the Dieterich (1994) rate-and-state model, which assumes a population of faults above steady-state (initially already on their way to failure), can be approximated with such a threshold Coulomb model due to the tendency of aσ to act as a stress threshold. On the other end, the application of the rate-and-state model to a population of faults below the steady-state regime also results in introducing a threshold in the rate-and-state model as well (Heimisson et al., 2022), accounting for the population of earthquake sources that are initially far from instability which is assumed negligible by Dieterich (1994). In this case, the question remains whether C cpt is indeed solely representing the initial stress state, or rather acting as a proxy variable that also encompasses effects of time-dependent nucleation. (Table 2) To address these questions, we consider a Coulomb model with a stress threshold representing the initial strength excess on the triggered faults. The Coulomb model is formulated as follows: where V is the representative volume over which seismicity is recorded, α c is a scaling factor defined as the change in pore pressure per slip event per unit volume (Nur & Booker, 1972), and f c is the probability density function representing the distribution of threshold triggering pressure needed for the Coulomb stress change to exceed the initial strength excess. Following the observation that poroelastic stress changes are minimal compared to pore pressure changes, they are ignored hereafter for simplicity. The derivation of Equation 16, which is the time derivative of Equation 7, is given in Appendix A of Segall and Lu (2015). The integral is restricted to where stress changes are positive, and to account for the Kaiser effect, the integral is further limited to where the past maximum pore pressure has been exceeded. Following Bourne et al. (2018) and Smith et al. (2020), we next assume a population of faults with randomly distributed strength excess using a formulation that has been found to provide a good model of seismicity induced by gas extraction from the Groningen gas field. Seismicity starts once the Coulomb stress change exceeds the lowest value of the initial strength distribution. According to the extreme value theory, the tail of the distribution can be represented by a Generalized Pareto distribution, leading to an exponential increase of seismicity for a constant loading rate (Bourne et al., 2018). This general formulation is valid to simulate the onset of seismicity but it does not allow the transition to a steady state regime where seismicity rate would be proportional to the loading rate. We therefore assume a Gaussian distribution of initial strength to allow for the transition to steady-state (Smith et al., 2020), and express it in term of the distribution of threshold pressure where θ 1 and θ 2 are the mean and standard deviation of the distribution, respectively. The best fitting model is found with respect to θ 1 and θ 2 within the range of 0.01-5 MPa for both parameters. α c is adjusted to match the total number of events, much like r b of the rate-and-state model. This model is hereafter referred to as the threshold Coulomb model.
The model fit in time and space are shown in Figures 6e and 8e, respectively, with θ 1 = 0.66 MPa, θ 2 = 0.28 MPa, and α c = 14.3 kPa/event ⋅ m 3 . The model fits the observations well in time, with a KS-statistic of 0.029 but significantly overestimates the extent of seismicity in space, which was also a main issue with the standard Coulomb failure mode ( Figure S7 in Supporting Information S1). The model is also less sensitive to rapid variations of the injection rate compared to the rate-and-state models, with relatively muted changes in the seismicity rate in-between injection cycles. Such sensitivity is seen to grow with the time scale of stressing rates; Figure 12 shows the response of the both the Coulomb and rate-and-state models with the duration of injections and pauses multiplied by factors of 0.1 and 10 (parameters are fixed to those that produced Figures 6d and 6e). While both models show more rapid variations of the seismicity rate relative to the injection rate for longer injection duration, the tendency is significantly greater in the threshold Coulomb model. For longer injection duration, the models show rather good agreement between each other although the threshold Coulomb model predicts lower t r with increasing time. Similar sensitivities may be observed with respect to the choice of θ 1 . While both the Coulomb and rate-and-state models may provide sufficient hindcasting tools for the same observation, the calibrated models produce very different forecasts for injection scenarios with duration of injection different from those used for calibration. In addition, they may produce different predictions in space for similar predictions in time. The comparisons suggest that the stress state with respect to failure and nucleation effects must be modeled separately, as done for example, in the threshold rate-and-state model of Heimisson et al. (2022), especially for fast injection cycles commonly employed in EGS operations where the effect of delayed nucleation may not be appropriately represented by the inclusion of a stress threshold in Coulomb models.
We remark that our modeling allows estimation of the best fitting values of a to between 0.00006 and 0.0002, which is significantly lower than the values inferred from laboratory measurements, generally ranging between 0.01 and 0.001 (Marone, 1998). Yet, the importance of rate-and-state effects in matching the observations in both space and time suggest that even such low values do not yield, for the injection schedule studied here, the rate-independent behavior that could be matched with a threshold Coulomb model. The reconciliation of field-inferred values of aσ and laboratory measurements is still paramount for eventual application of such models toward seismicity forecasting. One possible explanation is that spatial heterogeneities lead to elastic interactions that produce globally inferred values lower than that in a homogeneous equivalent (Dublanchet et al., 2013). It is also important to note that the model of Dieterich (1994) is a rather limited representation of the full complexity of rate-and-state friction. For example, the model simulates a population of spring-slider nucleation sources, whose qualitative differences in their behavior to more realistic finite fault models have been displayed for numerous aspects of rupture characteristics. Additionally, the model neglects the effect of variable effective normal stress on nucleation size, as the number of active nucleation sources remains constant throughout (Alghannam & Juanes, 2020). Further development of the model with a more holistic representation of rate-and-state friction would prove valuable for induced seismicity forecasting.

Origin of Omori-Law Decay Following Hydraulic Stimulation
The rate-and-state model reveals that the post shut-in Omori-law decay at Otaniemi depends strongly on the stress relaxation process by pore-fluid diffusion and cannot be explained solely by nucleation effects. The dependence on both nucleation and stress relaxation can be demonstrated by a sensitivity analysis of the relaxation timescale of the Omori law, t r , to parameters a, the rate-and-state friction parameter and k, the permeability. We find the most direct relationship to be that between the ratios of t r and the characteristic diffusion time, = 2 , to t a as shown in Figure 13 where t r is measured by fitting the Omori law to shut-ins following single boxcar injections under the rate-and-state model. Thus, t r is more strongly dependent on t c . The positive relationship t r and t c follows the intuitive reasoning that higher diffusivity would result in more rapid relaxation of induced pressure and consequently a faster decay of the seismicity rate. Our observations are consistent with the suggestion that the empirical Omori-law would be a result of stress relaxation by pore pressure diffusion (Almakari et al., 2019;Langenbruch & Shapiro, 2010;Miller, 2013). This explanation seems certainly reasonable in the context of EGS stimulations where pore pressure variations are particularly large.
The dependence on stress relaxation implies that t r also depends on injection duration ( Figure 13). The sensitivity analysis is performed with a and k fixed at 0.001 and 10 −16 m 2 , respectively, while the injection duration varies between factors of 0.1-100 of t c . The plot shows a nonlinear relationship between t r and the injection duration, t I , with an initial increase followed by a decrease. The trend exhibits a strong correlation with the seismicity rate at the time of shut-in. For shorter injections, the seismicity rate continuously increases prior to shut-in, increasing the time to relax to background levels. This is until the seismicity rate begins to decrease for continued injection, as pore pressure reaches steady-state conditions, and further nucleation is suppressed by the Kaiser effect ( Figure  S5 in Supporting Information S1). Consequently, t r decreases as well, as it takes less time to relax the lower seismicity rate. A similar effect arises due to the finiteness of the computational domain-the further distances where the seismicity rate would continue to increase at later times are cut-off. The sensitivity of t r to the total injected volume is consistent with the observation that the Omori law relaxation time at shut-in increases with time during the EGS stimulation at Otaniemi (Avouac et al., 2020).

Application of Models to Seismicity Forecasting
The models so far have only been applied in a hindcasting sense such that the data has been used in its entirety in order to tune the model parameters. We test the ability of the models to truly forecast induced seismicity in Otaniemi by limiting the range of the data used for training the models. Forecasts from the best fitting physical model (rate-and-state model with c true = c bu - Figures 6d and 8d) and the spatio-temporal convolution model are shown in Figures 14 and 15, respectively. The rate-and-state model is able to produce a forecast comparable to the hindcast using just the first injection stage as the training period with a similar value of a = 0.00005 although with significantly lower = 0.1 kPa/yr and r b = 0.39 events/day. With the same training period, the convolution model performs rather poorly, largely due to the estimation of t r at the end of first injection stage substantially lower (2.9 hr) than the average value throughout the entire injection schedule. The forecast is significantly improved by including the second injection stage within the training period, which now consists of the Omori decay observed during the injection pause at around the 450 hr mark that significantly increases the estimated value of t r to 10.4 hr.
It is likely that the rate-and-state model is more robust to the length of the training period than the convolution model due to the fact that c true is fixed at c bu which matches the pressure history of Otaniemi in its entirety ( Figure 7b). As discussed in Section 8.2, diffusivity plays a significantly stronger role in governing the rate of Omori decay than the tuning parameters of the rate-and-state model, namely a and . Thus, the rate-and-state model seems suited to perform well in forecasting applications given an accurate estimation of the diffusivity. Forecasts from the convolution model could also be improved by accounting for the increase in t r with cumulative injected volume as observed in Otaniemi (Avouac et al., 2020).

Influence of the Kaiser Effect
We have seen that the fit to the temporal evolution of seismicity is improved when the Kaiser effect is reset at each new stimulation stage. Although the clouds of seismicity generated during each stimulation stage overlap largely (Figure 3), this reset is justified as each new stage implied the stimulation of a new volume near the wellbore. Without such an adaptation, the seismicity rate is predicted to significantly lower during the second half of the injection history ( Figure 6c) along with large regions of seismic quiescence near the injection well (Figure 8c). This also implies that the efficacy of the convolution model-which does not account for the Kaiser effect at all-depends strongly on the apparent absence of the Kaiser effect in Otaniemi.
The physical mechanism behind the activation of new volumes is unclear given the diffuse and rather random structure of the relocated catalog (Figure 3). If this argument is dismissed based on relocation uncertainties, one could pose that a low diffusivity stimulated non-overlapping volumes from one stage to the other. However, such a low diffusivity should manifest in inconsistencies with the observed catalog in time, for instance, a longer apparent relaxation time during shut-ins. Rather, the need to reset the stressing history for the models to reproduce the observations in Otaniemi more likely implies the creation of new hydraulic pathways due to the fracturing Figure 13. Dependence of Omori law decay on fluid transport properties: t r of Omori law decay in response to single boxcar injections under the rate-and-state model are plotted in terms of t c and t a (left). t r , shows a stronger dependence on t c , or the diffusivity, than on t a . Namely, longer diffusion times result in longer relaxation times of the seismicity rate. t r also shows strong dependence on injection duration, t I (right). t r first increases with increasing seismicity rate at time of shut-in, before decreasing as steady-state stress conditions are reached when the seismicity rate decreases as well due to the Kaiser effect ( Figure S5 in Supporting Information S1).
nature of the stimulation that activated new nucleation sources (Cladouhos et al., 2016). Such phenomenon would depend on both the physical properties of the injected medium such as its fluid transport properties and fracture toughness, and the injection scenario, especially any spatial variation of the injection location.

Validity of the Convolution Model
Our study shows that, in the context of the Otaniemi injection schedule, the seismicity response to injections in time and space can be approximated with a simple convolution model. This model ignores all the sources of nonlinearity that may arise from the coupling between fluid flow and deformation, the earthquake nucleation process, the initial strength distribution and Kaiser effect. It is therefore not obvious that this approximation would be applicable to other induced seismicity context or for other injection schedules. We have therefore used our physical model to explore the parameter regimes under which the linear convolution method is able to match the rate-and-state model. The results are presented in Text S6 in Supporting Information S1. We found the success of the convolution model to depend strongly on the impact of the Kaiser effect on the linearity of stress evolution for the given injection schedule although it is also seen to be robust to nonlinear effects from delayed nucleation.

Conclusion
Physical models based on rate-and-state friction and stress changes due to pore-pressure diffusion and poroelasticity can successfully reproduce the seismicity observed during the EGS simulation which were carried out on the Otaniemi campus near Helsinki, Finland. While pore pressure measurements at the well indicate two possible diffusivities that fit either the build-up of pressure or its drawdown, the physical model suggests that the diffusivity of the medium is likely closer to the higher diffusivity fitting the build-up. We find that the effect of time-dependent nucleation is crucial in reconciling the higher diffusivity with the spatio-temporal distribution of triggered seismicity. Namely, the tendency of the parameter aσ to act proportionally to a triggering threshold significantly affects the apparent diffusivity inferred from the triggering front in Otaniemi. However, the effect of nucleation cannot be approximated well by introducing a stress threshold in the standard Coulomb friction model, at least in the context of rapid variations of injection rates common in EGS operations. We remark that there are significant portions of the relocated catalog that the models do not fully capture in space, such as the back-propagation front or far-field seismicity, although a significant portion of the observed far-field seismicity may have been due to leaks in the well casing. The Omori law decay observed in Otaniemi is seen to depend Figure 15. Partial forecasting of induced seismicity by convolution model: Ability of the convolution model to forecast induced seismicity is tested by limiting the portion of the data used for model tuning. The top two rows compare forecasts using the first one and two injection stages as training periods where t r is estimated to be 2.9 and 10.4 hr, respectively. The forecast using solely the first injection stage as the training period significantly underestimates t r and underpredicts the seismicity rate for the rest of the injection history. The forecast using the first two injection stages as the training period is comparable to the hindcast of Figures 6a and 8a, with only a marginally higher KS-statistic of 0.047 and lower log-likelihood of 175,430. strongly on fluid transport properties in the physical model. Lastly, the physical model indicates that the Kaiser effect is present in Otaniemi, weakened by the successive variation of injection locations between different stages.
We show that a statistical model whereby the seismicity rate is predicted in time and space by convolution of a kernel function inspired by Omori law decay with the injection rate can provide a good match to the seismicity observed in Otaniemi. The existence of such linear convolution kernels is consistent with the identification of systematic decay patterns in the rate densities calculated by adaptation of the cascading algorithm of Marsan and Lengline (2008) to induced seismicity. The statistical model is extended to space by incorporation of a half-norm distribution component to the kernel mimicking the behavior of the physical model. We find that the validity of the method, which assumes a linear relationship between the injection history and the induced seismicity rate, depends strongly on the presence of the Kaiser effect. The convolution model would be applicable toward injection schedules that minimize the impact of the Kaiser effect by decreasing injection durations relative to the local diffusion time or by variation of injection locations in space.
The physical model presented in this study makes a number of assumptions. One assumption is that it is appropriate to use Darcy's Law, which was established for a homogeneous porous medium, to model the flow in the fractured crystalline bedrock. Although the assumptions largely stem from the lack of data on local heterogeneities or anisotropy, neglecting presence of vertical or horizontal geological layers may be appropriate for Otaniemi where the objective is to fracture a largely crystalline medium. The model also ignores the effect of pore-pressure change on permeability. This is clearly an oversimplification as, in the case of fractured flow, the permeability increases substantially with pore pressure (Acosta & Violay, 2019;Cappa et al., 2006;Cornet & Jianmin, 1995;Evans et al., 2005). Common values of in-tact granite under comparable pressure are documented to be closer to 10 −21 m 2 (Brace et al., 1968), several orders of magnitude lower than that of the best fitting model (10 −16 m 2 ). Indeed, there are indications of changes in the diffusivity from the evolution of the injectivity index, or the ratio of injection rate to the well-head pressure ( Figure S10 in Supporting Information S1). Periods of heightened injectivity are well-correlated with periods of high seismicity rates, likely due to seismicity-induced increase in permeability. Reconciling the full scope of pressure variations at the well and the spatio-temporal patterns of observed seismicity would probably require an explicit account for the role of fractures and seismicity on permeability. Lastly, stress perturbations due to thermoelasticity can also be significant for EGS operations where temperature gradients between the injected fluid and surrounding medium are large (e.g., Gens et al., 2011;Im et al., 2017;Rutqvist & Oldenburg, 2008).
The modeling methods presented here could be useful in designing EGS operations or to interpret induced seismicity observations in terms of transport properties within the stimulated volume. They could additionally serve as a basis for a probabilistic TLS or be incorporated in a control and optimization framework such as the one presented by Stefanou (2019). At the 'moment, TLS are deterministic and based entirely on the observed maximum magnitude (Ader et al., 2020;Bommer et al., 2006;Kwiatek et al., 2019;Majer et al., 2007). As such, a red light event can be triggered by the occurrence of a rare event, with improbably large magnitude, that might not necessarily reflect an increased hazard level. In addition, such TLS do not provide a way to anticipate the response to possible mitigation strategies. This is important as many operations have been terminated as the original TLS design proved to be insufficient in preventing "red-light" incurring events (Grigoli et al., 2017;Majer et al., 2007). To alleviate that issue, our forecasting methods could for example, be incorporated in "Adaptive Traffic Light Systems" (Wiemer et al., 2015), which are based in a real-time assessment of probabilistic hazard.