A Coulomb Stress Response Model for Time‐Dependent Earthquake Forecasts

Seismicity models are probabilistic forecasts of earthquake rates to support seismic hazard assessment. Physics‐based models allow extrapolating previously unsampled parameter ranges and enable conclusions on underlying tectonic or human‐induced processes. The Coulomb Failure (CF) and the rate‐and‐state (RS) models are two widely used physics‐based seismicity models both assuming pre‐existing populations of faults responding to Coulomb stress changes. The CF model depends on the absolute Coulomb stress and assumes instantaneous triggering if stress exceeds a threshold, while the RS model only depends on stress changes. Both models can predict background earthquake rates and time‐dependent stress effects, but the RS model with its three independent parameters can additionally explain delayed aftershock triggering. This study introduces a modified CF model where the instantaneous triggering is replaced by a mean time‐to‐failure depending on the absolute stress value. For the specific choice of an exponential dependence on stress and a stationary initial seismicity rate, we show that the model leads to identical results as the RS model and reproduces the Omori‐Utsu relation for aftershock decays as well stress‐shadowing effects. Thus, both CF and RS models can be seen as special cases of the new model. However, the new stress response model can also account for subcritical initial stress conditions and alternative functions of the mean time‐to‐failure depending on the problem and fracture mode.

The CF model assumes instantaneous earthquake triggering when the absolute Coulomb stress exceeds a static strength threshold. The predicted earthquake rate correlates linearly with the stressing rate if a pre-existing population of potential faults with uniformly distributed stress is assumed, as long as stress increases. In contrast, a complete absence of seismicity is expected during periods of stress shadows after stress drops. Because a retarded triggering of earthquakes is not considered, the CF model cannot explain aftershocks solely by the coseismic stress changes of the mainshock. It requires the involvement of secondary stress changes to reproduce the characteristics of aftershock sequences, such as delayed stress changes due to (a) poroelastic response (Nur & Booker, 1972), (b) pore-pressure diffusion (Miller et al., 2004), or triggered postseismic slip, the so-called afterslip (Perfettini & Avouac, 2004).
The RS model posited by Dieterich (1994) is based on the assumption that fault slip is described by a specific constitutive RS dependent friction law derived from laboratory experiments. It assumes that frictional instabilities occur in a population of fault patches in a way that the earthquake rate is initially constant for constant tectonic stressing. In contrast to CF, it depends only on the stress change and not on the absolute stresses. RS explains, like CF, a constant earthquake rate for constant stressing and a rate decrease during periods of stress shadows. Additionally, it reproduces an Omori-Utsu type aftershock decay following a positive stress step. While the CF model has additional to the stress loading and threshold only one model parameter, the RSM has three model parameters. Its derivation is quite complex and uses several additional intrinsic assumptions. In a recent study, Heimisson and Segall (2018) revisited the RS theory with a different assessment and interpretation of the various model assumptions.
The new effective media model developed in this study is a modified CF model, which accounts for delayed nucleation of both tensile and shear cracks, simply by a stress-dependent mean time-to-failure. We first introduce the concept of the effective media approach and the modified CF model, derive the theoretical implications, and then discuss the model predictions in context with some field observations for induced seismicity, aftershock decays, and stress shadow effects in comparison to the established CF and RS models.

Concept
The earthquake rate is the number of events in a given volume and time interval with magnitudes above a completeness value M c . We define an elementary rock volume V that undergoes a uniform time-dependent change in stress or pressure. The volume is cut through by pre-existing faults and fractures, which may possibly support an earthquake rupture (Figure 1a for illustration). The faults or neighboring rock contains local stress peaks ("sources") of uneven areas of varying size, which may be viewed as asperities. In the effective medium approach, the interaction of the sources is considered by replacing the rock volume with a substitute medium with average properties (effective modules) depending on the fault density (Dahm & Becker, 1998). The distribution of average background stress (represented by principal effective stresses ′ 1 and ′ 3 ) and Coulomb stresses at the individual sources can be visualized in a Mohr Coulomb diagram (Figure 1b). Although the absolute stress at each source within V can vary and also be different from the ambient background stress, the change in Coulomb stress is assumed to be equal. For instance, the Coulomb stress changes Δσ may occur by a redistribution of stress in V after a major earthquake, by other types of internal aseismic dislocation sources beside V or by pore pressure changes in the rock volume itself. In the CF model, an earthquake is instantaneously triggered when the Coulomb stress σ c at the source exceeds the inherent cohesive strength S 0 (Figure 2). For a constant and uniform stressing rate , the time-to-failure t f for a source is in this case simply given by = ∕̇ with ζ = S 0 − σ c being the difference between strength threshold and the actual Coulomb stress at the source. Thus, a constant rate r 0 results in the CF model for a uniform pre-stress density distribution ( ) = ( 0∕̇) ( ) , where H is the Heaviside function (H(x) = 1 for x ≥ 0 and 0 else).
While the deterministic failure criterion of the CF model is simple and thus attractive, it is an oversimplification. Instead of a fixed threshold and instantaneous triggering, instabilities are realistically expected to occur at various stress levels with stress-dependent nucleation times for instabilities to grow into seismic ruptures. Based on experimental data on stress corrosion of rocks (Innocente et al., 2021) and subcritical crack growth, the mean time-to-failure ̄ is an exponential function according tō The constant t 0 is the mean delay time for a critically stressed source (ζ = 0) and can be assumed very small. The parameter δσ controls the increase of ̄ with ζ and is in our study denoted as skin parameter (see Figure 2b). In the limit of δσ → 0, the time-to-failure becomes ∞ for ζ > 0 (stable) and 0 for ζ < 0 (instantaneous failure); thus leading to the CF model. Conceptual sketch to derive an effective media approach for the distribution of faults or asperities sources in an effective media approach. (a) An infinite number of sources of different length and orientation are randomly distributed in a rock volume V. As some of them have experienced slip, and others not yet, the stress field in the rock volume is heterogeneous. (b) The Coulomb Failure line is plotted in a Mohr stress diagram, where colored circles indicate the state of stress on individual faults. The heterogeneous stress in V was replaced by a substitute media with an effective, homogeneous stress. Each source is considered independently in the center of the rock volume with its orientation conserved. The Coulomb stress of each source, σ c , aligns on the Coulomb stress axis. If the rock volume is externally loaded (dashed Mohr circles), the same stress change is assumed to apply to all sources in V and the σ c is uniformly shifted on the Coulomb stress axis.

of 19
The exponential law has been suggested for both quasi-static crack propagation of tensile (Aktinson, 1984) and shear cracks in brittle rock (Ohnaka, 2013) if ζ is either interpreted as a function of the stress intensity factor or the difference between strength and Coulomb stress, respectively. For instance, Scholz (1968) used the exponential form of 1∕̄ to explain the mechanism of creep in brittle rock. The same law has been used to explain crack growth in ceramics (Wiederhorn et al., 1980) or the number of acoustic emissions under creep (e.g., Ohnaka (1983) and references therein). We follow Scholz (1968) and Ohnaka (2013) and define ζ by S 0 − σ c . However, the concept can be applied equally well to tensile crack seismicity or other forms of the mean decay times.
Given the mean time-to-failure ̄ , the mean failure rate for N sources is given by ∕̄ . Thus, the distribution of sources at the different stress levels has to be considered to calculate the total event rate of the volume. This time-dependent density distribution of sources is defined by χ = χ(ζ, t) and the total rate at time t becomes where ̄( ) is given by Equation 2 and independent of time. In order to use this relation for rate forecasts, the evolution of χ(ζ, t) is needed. To solve it analytically or numerically, we use the following assumptions and simplifications: 1. Each source in V is loaded by the same stress and acts independently (effective media approach).
2. In addition to external loading, each source is subject to static fatigue described by Equation 2. 3. Each source fails only once during a simulation, that is, the stress drop ΔS is significantly larger than the loading stress during the simulation, and failed sources can be removed from the distribution of available sources. 4. For simplifying the following analytical calculations, we ignore that ζ is limited by the maximum value S 0 .

Analytic Solutions
The model forecasts can be analytically derived for some simplified stress scenarios. These solution can be all derived from the total time derivative of χ being equal to the negative event rate, namely The detailed derivations of the following solutions can be found in Appendix A. . In the TDSR model, the probability is a smooth function due to a non-zero skin parameter δσ (red line). (c) Normalized distribution of source stresses are defined by an initial susceptibility function χ(ζ, t 1 ). The initial distributions at time t 1 for δσ = 0 are assumed uniform at ζ > 0 (greyish area, CF model). For δσ > 0 the equilibrium distribution is a smooth function (bluish area). Due to stress loading or unloading, the distributions are shifted either to higher or lower levels on the stress axis with time.

Uniform Initial Distribution With Constant Loading
For the case of an initially uniform stress distribution for ζ > ζ min at t = 0, that is, χ(ζ, 0) = χ 0 H(ζ − ζ min ) (gray line in Figure 2c) and a constant stressing rate for times t > 0, the time-dependent seismicity rate is The solution can be compared to the rate expected in the CF and RS model. In the CF model, ( ) = 0̇( − min∕̇) , that is, the CF model, unlike ours, predicts a sudden onset of a constant seismicity rate at time = min∕̇ . In contrast, the original RS model of Dieterich (1994) always assumes a stationary background seismicity (r ≡ r 0 ) as starting condition and thus cannot deal with a subcritical initial stress state. Heimisson et al. (2022) extended the RS-framework for subcritical initial stress states (RS subcrit ). The solution of their Equation 1 for ΔS c = ζ min and a constant stressing rate yields ( ) = 0 ( − min∕̇) , that is, an instantaneous onset of a constant rate at = min∕̇ as the CF model.

Stationary Seismicity
The stress distribution associated to a constant seismicity rate R(t) = r 0 for a given stressing rate is given by (blue line in Figure 2c) Here, the subscript in χ s is used to denote that this distribution is stationary, that is, time-independent. The prefactor 0 = 0∕̇ expresses the susceptibility of the rock volume.

Stress Step
Here it is assumed that the initial pre-stress distribution is given by χ s related to constant stressing and seismicity rate, and r 0 , respectively. A stress step Δσ c is applied at time t = 0 followed by a loading with constant stressing rate at t > 0. The seismicity rate is in this case given by , and = 0̇ . While the CF model cannot predict Omori-type aftershocks, both equations are identical to the solutions of the RS-model for the same conditions, namely, Eq. (12) and (13) of Dieterich (1994) with Aσ ≡ δσ, r ≡ r 0 , ≡ , and ≡ .

Changing Stressing Rates
The solution for a changing stressing rate from for t ≤ 0 to the new constant stressing rate for t > 0 can be directly derived from Equation 7 by setting Δσ c = 0, yielding A ramp-like excitation in stress, or a box-car stressing rate with duration t b and for t < 0 and t > t b and elsewhere, can be constructed from Equation 9 by 6 of 19 where = −̇ln ( 1 − −̇) and R 1 given by Equation 9 and R 2 by Equation 9 with interchanging and .

Applications
While the analytic solutions (Equations 5-10) can be used for simple stressing histories, the integral in Equation 3 must be solved numerically for more complex cases. The numerical model implementation is straightforward, and the structure of the simple algorithm, which we used for the simulations presented in this work, is provided in Appendix A2.

Synthetic Examples
The following synthetic case studies are used to illustrate the response to some generic setups. Figure 3a shows the TDSR forecasts for the case of a uniform pre-stress distribution and constant loading . The numerical simulations of the TDSR model exactly resemble Equation 5. The three analyzed ζ min values of 0, 3, and 6 lead to different responses. For ζ min = 0, the rate decays with time until it reaches the stationary seismicity rate ∞ = 0̇ . In contrast, the seismicity rate is initially much lower than r ∞ and only starts to accelerate and converge with some delay for ζ min > 0. In comparison to the predictions of CF and RS subcrit , which are both identical for this setup, the onset of the seismicity is smooth and occurs earlier. The accelerated seismicity already reached the steady-state when the other models predicted the sudden onset. Thus, for a fixed onset time of the seismicity, the estimated ζ min value for TDSR will be larger than for CF and RS subcrit (e.g., see the application to Groningen in Section 3.2.1).

Aftershock Triggering
The TDSR model with stationary pre-stress distribution (Equation 6) leads to Omori-type aftershock sequences with p = 1 for positive stress steps at t = 0. Figure 3b shows the immediate rate increase in a double logarithmic plot as a function of time after the stress step, where the stressing rate remains constant, that is, = . The numerical simulations exactly match Equation 7. It is important to note that (a) the aftershock rate is directly proportional to the background rate r 0 , (b) the maximum seismicity rate at t = 0 + is r 0 exp(Δσ c /δσ), (c) the duration of the aftershock decay is given by ≡ ∕̇ , and (d) the time delay until the onset of the 1/t decay (c-value in Equation 1) depends on the stress step. The larger the step, the shorter the delay.

Stress Shadowing
The seismicity rate is commonly assumed to decrease during periods in which the absolute stress dropped below its previous maximum value (so-called stress shadows). This stress shadow effect is also known as Kaiser effect. To illustrate the response of the TDSR model to stress shadows, we run simulations with stationary pre-stress distributions. Figure 4a shows the result for sudden stress drops (Δσ c < 0), which are also described by Equation 7. For instance, a volume of rock on the rupture plane of a large earthquake may experience a co-seismic negative stress step when the rupture front has passed by. The TDSR model predicts identical to the RS model that the seismicity rate drops instantaneously in this case from r 0 to r 0 exp(Δσ c /δσ) and then recovers slowly until it reaches the background level at approximately 6 t a . Thus, the memory for stress drops is a few times larger than for positive stress steps, that is, aftershock sequences. In contrast, the CF model (i.e., the TDSR model in the limit for δσ → 0) predicts an absolute quiescence of length t a |Δσ c |/δσ followed by an instantaneous recovery of the background rate. Similarly, for cyclic loading shown in Figure 4b, the TDSR model predicts smooth transitions from almost zero rates to high activity and vice versa, while CF shows sharp transitions.

Arbitrary Stressing Histories in Space and Time
In general, stress changes are complex in time, for example, due to multiple stress sources such as earthquakes, aseismic slip, or pore-pressure changes, and no analytical solution is known. In addition, the stress evolution is usually spatially inhomogeneous because distances to sources vary and physical properties such as fault density and orientation or stress state are spatially heterogeneous (Hardebeck, 2021). In this case, the model must be numerically solved by discretizing the stress evolution in space and time and applying the algorithm described in Appendix A2 independently for each spatial grid point.
It is important to note that TDSR leads to identical forecasts for any arbitrary stress history as RS if the initial condition is critical, a reasonable assumption for natural seismicity at plate boundaries if no other information is available. This congruency is expected given that (a) the response to stress steps (Equations 7 and 8) and changes in the stressing rate (Equation 9) is for both models identical, and (b) any stressing history can be well In all time-dependent stress response simulations, an initially constant seismicity rate and a constant background stressing rate is assumed. While the rate-and-state model predicts the same response in all cases, the Coulomb Failure model (dashed lines) predicts a total quiescence of length t a |Δσ c |/δσ after the stress drop in (a and b) and during stress shadows (d) marked by horizontal lines in (c).

of 19
approximated by stress steps and periods of constant stressing rates using sufficiently small time steps. However, for confirmation, we also performed simulations with complex stress histories consisting of multiple stress steps with variable stressing rates in between, finding identical forecasts of TDSR and RS.

Applications to Real Data
Analytic solutions and simple synthetic model runs are helpful to clarify and demonstrate general model properties, but the ultimate goal of any seismicity model is the application to real data. In the following, we briefly discuss some examples of TDSR applications to induced and natural seismicity, which are related to the synthetic case studies.

Groningen
The Groningen field, Netherlands, is one of the largest gas reservoirs in the world. The gas production started in 1963, but the first earthquake related to gas extraction was only detected in 1991. At that time, the reservoir pore pressure dropped already approximately 10 MPa, related to a rock stress increase of similar magnitude, generally taken as evidence that before production, most Groningen faults were far from critical stress (Candela et al., 2018). Between 1991 and 2021, more than 350 earthquakes with magnitude M L ≥ 1.5 have been recorded.
The occurrence of felt earthquakes accompanying the gas production initiated dense instrumentation, numerous field studies, and the development and application of many seismicity models for explaining and quantifying the observed seismicity, see the recent review of Kühn et al. (2022). For example, a CF model version was developed by Bourne and Oates (2017) and Bourne et al. (2018) assuming that the stresses at the activated faults are in the tail of the pre-stress distribution. The latter assumption was used to apply extreme value statistics to explain the observed non-linear response to the stress changes. Furthermore, different RS model implementations have been applied to Groningen, among others Dieterich's model (RS) by Richter et al. (2020) and the subcritical RS model (RS subcrit ) by Heimisson et al. (2022). Recall that TDSR forecasts are identical to RS forecasts under the assumption of critical initial conditions. Thus, all Richter et al. (2020) results for RS can be interpreted as TDSR predictions with δσ = Aσ, r 0 ≡ r 0 , and ≡ ∕ using Equation 6 as initial condition. Previous model applications were based on estimated stress changes in space and time, including available fault information. Here, we do not want to compete with those detailed studies but only want to demonstrate the potentials of the TDSR model to account for subcritical stress states and leave a detailed analysis to future work. For simplification, we only consider the annual changes of the mean pore pressure in the field and assume that Coulomb stress changes are proportional to the size of the pore pressure drop.
We compare the TDSR forecasts with the corresponding results of the CF, RS, and RS subcrit models. To reduce the number of free model parameters, we constrained the pre-factors by the condition that the total number of forecasted earthquakes equals the observed events in 1991-2021. With this constraint, CF considering continuous tectonic stressing has no free parameter, CF assuming a subcritical initial stress state (CF subcrit ) has one (ζ min ), RS has two (Aσ and ), and RS subcrit has three parameters ( , and ζ min ). For both models, CF subcrit and RS subcrit , the stress gap had to be equalized by the gas production at the time when the seismic activity started; thus, we set ζ min = 8.5 MPa, which was equalized by the induced stress increase just before the first observed earthquakes in 1991. Furthermore, we set Aσ = 1 MPa which is in the range of the results of Richter et al. (2020). The free parameters of the TDSR model depend on the assumed initial conditions. Assuming a constant background seismicity, that is, a stationary TDSR pre-stress distribution according to Equation 6, leads to forecasts which are independent of t 0 and identical to the RS model for δσ = Aσ using the same stressing rate . However, the TDSR model also works for subcritical stress conditions with zero background stressing. Instead of , the model depends on t 0 and the parameters required to characterize the pre-stress condition; for example, ζ min for a uniform pre-stress distribution, or the mean ̄ and standard deviation ζ σ for a Gaussian distribution. Figure 5b shows the results of the TDSR model in comparison to the other models, where we set the remaining free parameters in order to fit the observed annual rates of earthquakes with magnitudes larger than 1.45, that is, the estimated completeness magnitude since 1993 (Dost et al., 2017;Heimisson et al., 2022). We find that 1. The CF model with uniform pre-stress conditions cannot fit the observations, neither for critical nor for subcritical conditions. 9 of 19 2. TDSR and RS assuming an initially constant background seismicity rate lead to identical fits. For a reasonable data fit, the background stressing rate must be extremely low ( = 3.3 Pa/year in Figure 5). 3. RS subcrit leads to almost the same result as RS but with a significantly larger background stressing rate, namely 17 kPa/year for the same Aσ value. 4. To get similar fits, the pre-stress conditions in the TDSR model depend on t 0 ; the shorter t 0 , the larger the pre-stress gaps, or vice versa. For a mean nucleation time at critical stress of t 0 = 10 −4 years (≈53 min), the stress gap is significantly larger (ζ min = 22 MPa, respectively ̄ = 24.5 MPa) compared to 8.5 MPa estimated by CF subcrit and RS subcrit . 5. The shape of the ongoing seismicity rate decay depends strongly on the assumed pre-stress distribution shown in Figure 5a. While the predictions of the RS models and the TDSR model either with background stressing or uniform pre-stress distribution behave similarly, the decay is much steeper if a Gaussian distribution of sources is assumed. In the latter case, the depletion of available sources starts to dominate the activity at presence and is expected to lead to a rapid disappearance of the activity in the near future. Because this could have important implications for seismic hazard studies, a more detailed analysis is needed using more specific stress calculations in space and time and an estimate of the initial stress distribution considering the regional stress field and fault orientations. Such a comprehensive analysis is beyond the scope of this study and is the subject of future work.

Fits to Aftershock Sequences
Important applications of seismicity models are forecasts of aftershock sequences because aftershocks are sometimes very destructive and deadly. Major aftershock sequences are typically observed at tectonic plate boundaries, where a constant background stressing rate can be assumed for loading the fault system. Thus, in the absence of The resulting seismicity response of the TDSR model to the mean pressure changes in the Groningen gas field, using t 0 = 10 −4 years and δσ = 1 MPa. Here, points refer to the annual rate of observed earthquakes with magnitudes larger than 1.45, that is, the estimated completeness magnitude since 1993 (Dost et al., 2017;Heimisson et al., 2022). The largest earthquake occurred 2012 with a magnitude 3.6. For comparison, the corresponding curves are shown for the subcritical Coulomb Failure (CF) model (CF subcrit ) and the subcritical rate-and-state (RS) model of Heimisson et al. (2022) (RS subcrit ). In both cases, an initial stress gap was set to 8.5 MPa, reached just before the first observed earthquakes in 1991. In addition, a background stressing rate of = 17 kPa/year is assumed in RS subcrit . Note that the RS model of Dieterich (1994) leads to an identical forecast as TDSR s .
10 of 19 any additional information, a constant background rate as initial condition is reasonable in this case. Then, as discussed before, the TDSR model forecasts are identical to the RS model with δσ = Aσ.
Thus, the detailed existing studies and comparative tests for the RS model can be exploited in terms of the TDSR model. Those model studies made use of detailed slip models for the mainshocks and large secondary stress changes to estimate the seismicity rate changes in space and time; among others, for the 1992 M7.2 Landers sequence (Hainzl et al., 2009), 2004 M6 Parkfield and the 2011 M9 Tohoku sequences (Cattania et al., 2014), the 2010-2012 Canterbury sequence (Cattania et al., 2018), and the 2016-2017 Central Italy earthquake cascade (Mancini et al., 2019). Those model applications provided, in general, good data fits. Furthermore, pseudo-prospective forecast experiments were conducted for several sequences using the consistency and comparative tests introduced by CSEP (Collaboratory for the Study of Earthquake Predictability, https://cseptesting.org). The comparison of the TDSR/RS forecasts with those of the purely empirical ETAS model showed that the physics-based models partly outperform the empirical ones (Cattania et al., 2018;Mancini et al., 2019;Woessner et al., 2011). These studies for aftershock sequences also provide estimates of Aσ, that is, δσ, in the range of 0.01-0.1 MPa.
As mentioned in Sections 2.2.3 and 3.1.2, TDSR predicts that the duration t a of the aftershock decay is given by = ∕̇ . Thus, it is inversely related to the tectonic (background) stressing rate , implying that aftershock sequences in intraplate regions last longer. For a stress drop of 1 MPa, assuming a reloading time of 100 and 1,000 years corresponds to stressing rates of 0.01 and 0.001 MPa/year, respectively, which would result in t a in the range of 1-10 years and 10-100 years, respectively, using the estimated range of δσ values. These values are consistent with empirical observations, for example, the results of Hainzl and Christophersen (2017) for empirical fits to more than 300 recorded earthquake sequences, which yielded a mean value of about 30 years, but with a large proportion of significantly shorter t a -values.
Another prediction of the TDSR model is that the time delay until the onset of the 1/t decay (c value) depends on the size of the stress step. The larger the step, the shorter the delay, according to c = t a exp(−Δσ c /δσ). Therefore, the onset of aftershock decay near the rupture should be quicker than for aftershocks in the far-field. For example, assume t a in the range of 1-100 years and the mean value δσ = 0.05 MPa. Then, a stress change of 1 MPa results in c values between 0.06 and 6 s, while a stress change of 0.1 MPa results in c values between 50 days and 13.5 years. This example shows that the predicted difference is enormous and could have far reaching implications. In practice, however, the test is difficult. First, short c values cannot be detected due to the incompleteness of earthquake catalogs in the first hours to days after a mainshock when small earthquakes cannot be recorded because of overlapping waveforms (Hainzl, 2016). Second, longer c values are masked by secondary aftershock triggering. More distant aftershocks are often triggered by closer, early aftershocks rather than the more distant mainshock. Therefore, analyzing and testing the predicted stress dependence of the c value requires careful consideration.
Finally, TDSR predicts that the response to stress steps is proportional to the background rate r 0 (Equations 7 and 8 and Figure 3b). This prediction is in contrast to the empirical ETAS model, where the aftershock productivity is assumed to be independent of the background rate. Thus, TDSR predicts no aftershocks in ductile regions deforming aseismically and low aftershock productivity in regions of low seismic coupling. The latter agrees with the observation of Hainzl et al. (2019) that aftershock productivity is proportional to the seismic coupling derived from geodetic inversions in the Northern Chile subduction zone.
As seen by Figures 3b and 4a, TDSR predicts an Omori decay for positive stress steps and quiescence for stress drops. A combination of both is predicted for stress distributions with a negative mean value but a tail including positive values. Helmstetter and Shaw (2006) and Marsan (2006) theoretically showed that the aftershock decay in this case transitions to a quiescent phase after a short time. Such a distribution of stress steps is expected near the mainshock rupture, where stress drops on average but small-scale heterogeneities with positive stress changes exist. The prediction is consistent with observations for megathrust earthquakes, as recently shown by Toda and Stein (2022).
In general, mainshocks may also affect subcritically stressed regions. TDSR forecasts assuming subcritical prestress conditions differ from forecasts of the standard RS model and its subcritical version. While the response depends on the details of the prestress distribution, we explored a simplified setup to analyze the general effect. We compared the response to a stress step of Δσ c /δσ = 10 assuming critical initial conditions (uppermost curve 11 of 19 in Figure 3b) with the aftershock activity for subcritical prestress distributions where the initial values of χ s in Equation 6 are set to zero for ζ < ζ min . We find that values of ζ min /δσ between 0 and 8 resemble the curves in Figure 3b related to reduced stress steps Δσ c /δσ varying between 8 and 2. Thus, for increased subcriticality (ζ min / δσ), the aftershock productivity reduces, and the c-value increases. This effect could explain the low aftershock productivity (high Båth's law value) observed in stable intraplate regions (Page et al., 2016), where subcritical stress conditions are generally expected.

Examples of Stress Shadows Effects
Induced earthquakes are well suited as case studies for stress shadows when they are associated with a controlled, short-term variations in Coulomb stress, for instance related to periodic fluid injections in basement rocks or cyclic thermal loading of a mine gallery.

KTB Experiment
The first example treats a sequential fluid production and injection experiment in 4 km depth at the German Continental Deep Drilling site (KTB) in Windischeschenbach, Bavaria, Germany ( Figure 6). The KTB consists of two accessible boreholes drilled 200 m from each other into crystalline rock, the 4,000 m deep pilot well and the 9,100 m deep main well . Two distinct fault systems cut through the hydraulically connected wells. The hydraulic experiment started in June 2002 with a one-year production of 22,300 m 3 of saline water of 120°C from the open section of the pilot well between 3,850 m and 4,000 m depth . The drawdown of the fluid level in the pilot well reached 605 m at the end of the production phase in June 2003. The fluid level in the main well gradually dropped from zero to 50 m below the surface during the production phase, indicating a hydraulic connection between the two wells likely originating at a leakage in the casing of the main well between 5,200 and 5,600 m depth (Grässle et al., 2006). In June 2004, after a one-year recovery phase following the production, an injection experiment was conducted until April 2005. Overall, 84,500 m 3 of freshwater, about four times the produced volume, was re-injected during 10 months in the same open-hole section

of 19
of the pilot hole. Although the injection rate was more or less constant between 185 and 196 l/min, the wellhead pressure varied between 9 and 12 MPa . In October 2004, 114 days after the start of the injection, the main hole became artesian with a rate of ∼6.9 × 10 −4 m 3 /min (∼1 m 3 /day) . Because of the hydraulic connection between the two wells, the water level in the main well directly measured the time-dependent pressure change ΔP f at the leakage point at a depth of ∼5,500 m. Using this and other information on the local hydraulic structure, Grässle et al. (2006) modeled the pressure field between the two wells using a 3D diffusivity model. As we are interested in pressure changes only, our modeling approach is simpler and involves a poroelastic diffusion modeling in full space (Rudnicki, 1986). A hydraulic diffusivity of D = 0.033 m 2 /s is found to explain the ΔP f variations at the open section of the main well (Figure 6d). Here we used a rigidity of 1 GPa, a drained and undrained Poisson ratio of 0.25 and 0.3, respectively, and a Biot-Willies constant of α = 0.1. The model is used to extrapolate into the artesian phase during injection, where no water level measurements were possible.
Induced seismicity was recorded at a borehole geophone in the main hole (three components, 15 Hz, 1,000 Hz sampling rate, clamped at different depths between 1,950 and 3,000 m) and a surface array of five short-period stations in a distance of up to 3 km from the main hole (Haney et al., 2011;Shapiro et al., 2006). The largest 104 induced earthquakes had magnitudes between −1.8 ≤ M L ≤ 1 and could be accurately located from both networks. The events occurred between 3.3 and 4.9 km depth near the main well (Figures 6a and 6c). The borehole geophone was very close to the center of the cluster. In total, about 3,040 event arrivals were detected from which 2404 clustered micro-earthquakes were located with magnitudes M L ≥ −3.8 and with a magnitude of completeness of M c ∼ −2.75 (Haney et al., 2011). We use a more conservative limit of M c = −2.3 (Figure 6e).
The first seismicity occurred when the pore pressure at depth in the main hole reached and exceeded the original level from before 2002 (Figure 6d). We performed simulations to verify whether TDSR can explain the onset and time evolution of the earthquake rate. We assume that μΔP f dominates the Coulomb stress change, where μ is the friction coefficient. Our estimated Coulomb stress change thus covers the full experiment from the pumping over the recovery to the injection phase over altogether 1,200 days. As the experiment took so long, it poses a unique opportunity to test stress shadow models.
The stress state at KTB is known to be critical, where small stress perturbations triggers earthquakes (Zoback & Harjes, 1997). Thus, we assume a stationary pre-stress distribution according to Equation 6. The majority of the activity occurred with a distance between 300 and 600 m to the injection point, thus we averaged the TDSR-response function in this distance range and calibrated the resulting curve by the total number of observed earthquakes with magnitude exceeding M c . Figure 6f shows the result earthquake rates for δσ/μ = 0.9, 1.0, and 1.1 MPa using a tectonic stressing rate of = 30 Pa/year. Assuming μ ≈ 0.5, we estimate slightly larger δσ-values as for aftershock sequences. The onset and non-linear increase of the observed seismicity rate are well explained.

Morsleben
The second example involves the controlled refilling of an about 80 years old, abandoned salt gallery in the Morsleben mine, Germany, in a depth of 300 m below the surface. The salt-cement was backfilled continuously for approximately 182 days, with the exception of weekends and holidays, resulting in periodic variations of thermally induced Coulomb stresses in the roof region of the gallery. The periodic stress variations induced microseismic activity, so-called acoustic emissions in the frequency range of 1-20 kHz. A detailed description of the mine, the salt rock, the seismic monitoring system, and the location and magnitude estimation can be found in Köhler et al. (2009) andBecker et al. (2010). The high-quality catalog of AE events comprises hundreds of thousands of high-quality events which occurred in a confined volume above the backfilled gallery. Coulomb stress changes in the roof region were calculated with a 2D thermo-elastic finite element method using temperature measurements as boundary conditions (Becker et al., 2010). The observed AEs showed a high correlation to positive Coulomb stress changes but also clear stress shadow effects (Becker et al., 2010). Figures 7a and 7b shows a map-view and cross-section of the mine structure together with the seismic network and the AEs. Figure 7c shows the Δσ c changes in the selected rock volume of V = 15 × 15 × 10 m ∼ 2,250 m 3 in the southern region as defined in Becker et al. (2010). In this volume, Δσ c changes almost in phase with little spatial variations.
We applied the TDSR model using the modeled σ c (t). In contrast to the KTB case, the pre-stress is expected to be subcritical in the salt mine and tectonic stressing is zero. Thus, we assume a subcritical, uniform pre-stress distribution at the starting time, that is, χ(ζ) = χ 0 H(ζ − ζ min ). The susceptibility χ 0 is determined by the condition 13 of 19 that the total number of the predicted events should equal the observed ones. The remaining model parameters are ζ min , t 0 , and δσ. For simplicity, we set t 0 = 0.01 days but noticed that similar results are obtained for other t 0 values if ζ min is rescaled. Thus, the two parameters ζ min and δσ determine the shape of the predicted model rate. Figure 7d shows the comparison of the observed rates with the model forecasts for three skin parameters, δσ = 0.3 MPa, 0.4 MPa, and 0.5 MPa, and corresponding ζ min -values. TDSR shows a good fit to the data. The onset and width of the peaks as well as the relative heights of the cycles are well reproduced in all three cases. However, further increasing or lowering δσ leads to underestimating or overestimating the rate response to σ c (t) cycles. Thus, we can conclude that δσ is in the range between 0.3 and 0.5 MPa. This is almost the same range as for KTB, while approximately one order larger than the values estimated for aftershock sequences at active tectonic faults.

Discussion
Depending on the objective and application, various seismicity models are presently used, including statistical and physics-based models or combinations of both. A recent review of the strength and weaknesses of different models applied to the specific case of the Groningen gas field is given in Kühn et al. (2022). The most popular physics-based seismicity models are the variants of the CF and the RS models, while the ETAS model is the standard statistical model for describing short-term earthquake clustering.
The newly introduced TDSR model can be seen as a generalization of both physics-based CF and RS models. Specifically, the CF model is the limit of TDSR for δσ → 0, as discussed in Section 2.1. Thus, CF is the special 14 of 19 case of the TDSR for a vanishing skin parameter. Although CF can explain the occurrence of stress shadows (Kaiser effect), its prediction of sudden onsets and ends of total quiescences is not realistic. Furthermore, CF cannot explain Omori-type aftershock activity following sudden stress steps. This shortcomings vanish in the TDSR model with δσ > 0.
The relation of the TDSR model to RS is more complicated than to CF. As we have shown in this paper, the analytic and numerical results of RS are identical to the results of the TDSR model for the special case that the initial stress is in steady-state related to background stressing. In particular, the original RS model, as formulated by Dieterich (1994) and Heimisson and Segall (2018), uses the pre-condition = const > 0 and was not developed to study scenarios under = 0 . The assumption of initial steady-state conditions leads to a constrained pre-stress distribution. Heimisson et al. (2022) recently extended the RS approach to allow for subcritical starting conditions but still assumes non-zero background stressing rates, manifested in the model parameter = ∕̇ . In contrast, TDSR can be applied to arbitrary initial pre-stress conditions and does not generally require a constant background stressing rate. Simulations with ≈ 0 can be of particular interest for anthropogenic seismicity occurring in intraplate regions as in the case of our example in Morsleben; or, likewise, after driving a tunnel into a rock mass not subject to tectonic strain and stress rates.
While the predicted seismicity rates are identical for initial steady-state conditions related to > 0 , the concepts of both TDSR and RS are very different. The RS model assumes an infinite population of nucleation sites, where slip on each patch is described by a rate-and state-dependent constitutive law derived from friction experiments in the laboratory (Dieterich, 2007). For each of the isolated and non-interacting patches, a simplified spring slider model under friction is used to characterize the nucleation process by an interval of self-driven, accelerated slip, according to the aging law suggested by Ruina (1983). An earthquake is assumed to occur when the slip rate becomes very large or infinity. Finally, the seismicity model only depends on stress changes and not on absolute stress. In contrast, the TDSR model has a direct link to the absolute stress level by assuming that mean time-to-failure ̄ is a function of the absolute stress. We use an exponential function (Equation 2) that has been widely used to study subcritical crack growth and brittle failure in geological materials, both under tensional stresses (see review in Aktinson (1984)) and compressional (shear) stresses (Ohnaka, 1983;Scholz, 1968). Thus, it seems enigmatic why RS and TDSR result in the same solutions (for > 0 ), despite the contrasting concepts. The underlying reason might be that the RS approach also leads intrinsically to an exponential dependence of the time-to-failure on the absolute stress level. The derivation shown in Appendix B also provides the relation between t 0 and the microscopic and constitutive RS parameters. However, in general, the TDSR approach may also be used with other functional dependencies of ̄ on stress, for example, to study the occurrence of tensile cracks.
The simplicity of the TDSR model might offer additional possibilities. By explicitly tracking the stress distribution of sources, possible dependencies of the frequency-magnitude distribution on the stress state can be directly implemented. Particularly, the Gutenberg-Richter b-value is expected to depend on the absolute stress, with lower values (increased average magnitudes) for higher stress levels (Scholz, 2015). Such relations can be simply implemented and tested in the TDSR approach.
In contrast to the discussed physics-based models, the ETAS model does not rely on stress or stress changes at all. It is a statistical approach to describe short-term clustering of seismicity, particularly aftershock occurrence. For this purpose, it considers a stationary background rate and uses empirical relations to describe clustering: (a) the temporal aftershock decay (Omori-Utsu law), (b) the observed exponential dependence of the aftershock productivity on the mainshock magnitude, and (c) the spatial decay of the aftershock density with distance to the mainshock. Although, the physics-based model also derives all three properties (shown here and in Hainzl et al. (2010)) there are important differences: First, the background rate is additive in ETAS, while it is multiplicative in the physics-based models. It means that the short-term aftershock rate depends only on the mainshock magnitude in the ETAS model, while it also scales with the background rate in TDSR and RS. Second, the aftershock decay is infinite in the ETAS model, while it has a characteristic duration = ∕̇ in TDSR, which is better fitting empirical data on average (Hainzl & Christophersen, 2017). Third, the c-value of the Omori-Utsu relation is not, as in ETAS, a constant but depends on the stress step and is thus space-dependent in the TDSR model. Finally, and maybe most importantly, ETAS only accounts for activation and cannot explain stress shadows. The minimum rate at any location is the tectonic background rate. Despite these unphysical characteristics, the ETAS model mostly outperformed RS implementations in retrospective tests for major aftershock sequences 15 of 19 (Cattania et al., 2018;Woessner et al., 2011). The reason is likely that the physics-based models suffer from the large uncertainties of the mainshock-induced stress changes (Hainzl et al., 2009) and spatial heterogeneities of the physical properties, such as fault density, stress state, fault strength, and fluid pressure, leading to spatial earthquake clustering on small scales (Hardebeck, 2021). Additionally, ETAS accounts for secondary aftershock triggering by smaller magnitude events. However, in the case of anthropogenic seismicity, earthquake-earthquake triggering is less important, and seismicity is largely driven by the time-dependent stress changes related to the human activity. In the latter case, the physics-based models are superior, and ETAS might only be used to explain additional earthquake-earthquake triggering (Kühn et al., 2022).

Conclusions
Assuming simply an exponential dependence of the mean time-to-failure on the absolute stress, the proposed TDSR model can explain the most important characteristics of seismicity, namely aftershock triggering and stress shadowing. It is a generalization of the deterministic CF model where earthquakes nucleate instantaneously when the strength threshold is exceeded. Furthermore, TDSR leads to the same analytic and numerical solutions as the well-known RS model for the special case of critical initial pre-stress conditions related to a constant background rate, as can be assumed for natural seismicity in general. RS is based on a completely different concept, exploiting a rate-and state-dependent constitutive law derived from laboratory friction experiments. However, TDSR is not limited to initial steady-state conditions and can also simulate seismicity with subcritical initial conditions and zero tectonic stressing. The latter might be of particular interest for applications to human-induced seismicity.

A1. Derivation of the Analytic Solutions
Based on Equation 4 the following cases can be analytically solved.

A1.1. Time Evolution for an Initially Uniform χ in [ζ min , ∞]
For the case that the initial state is uniformly distribution according to χ(ζ) = χ 0 H(ζ − ζ min ) and the volume is loaded with constant stressing rate , the solution of Equation 4 is given by Inserting this solution in Equation 3 and solving the integral leads to

A1.2. Constant Rate for Constant Loading
Here we consider the case that the rate is constant, that is, R(t) = r 0 , for given stressing rate . In this case, the χ distribution is stationary, that is, d The constant c is determined by the condition ( ) = ∫ ( )∕̄( ) d = 0 , which yields = 0∕̇ considering that the integral of exp (− − − ) is equal to exp (− − ) ∕ , which is 1/a for x → ∞ and 0 for x → −∞ given a > 0.

A1.3. Response to a Stress Step Constant Loading
Here we consider a stress step Δσ c at time t = 0 followed by a constant loading rate for a source population which was initially at steady state related to a background stressing rate . In this case d d = −̇ and Equation 4 is equal to Given the initial condition χ(ζ, 0) = χ s (ζ + Δσ c ), the solution of this equation is given by Inserting this solution in Equation 3 and solving the integral leads to where it is again considered that the integral of exp (− − − ) is equal to exp (− − ) ∕ , which is 1/a for x → ∞ and 0 for x → −∞ given a > 0.

A1.4. Response to a Stress Step Without Subsequent Loading
For the specific case of = 0 , the solution of Equation 8 can be achieved by expanding the exponential term exp (−̇∕ ) of Equation A7 in a Taylor series and consider the limit → 0 .

A1.5. Change of Stressing Rate
For the case that at t = 0 the stressing rate changes ≠ , the solution provided in Equation 9 follows from setting Δσ c = 0 in Equation A7.