Seismic Diffusivity and the Influence of Heterogeneity on Injection‐Induced Seismicity

The spatiotemporal patterns of injection‐induced seismicity (IIS) are commonly interpreted with the concept of a triggering front, which propagates in a diffusion‐like manner with an associated diffusivity parameter. Here, we refer to this diffusivity as the “seismic diffusivity.” Several previous studies implicitly assume that seismic diffusivity is equivalent to the effective hydraulic diffusivity of the subsurface, which describes the behavior of the mean pressure field in heterogeneous porous media. Seismicity‐based approaches for hydraulic characterization or simulations of IIS using domains of homogeneous equivalent porous media are implicitly based on this assumed equivalence. However, seismicity is expected to propagate with the threshold triggering pressure, and thus not be controlled by the evolution of the mean pressure field. We present numerical simulations of fluid injection to compare the seismic and effective hydraulic diffusivities in heterogeneous formations (including fractured rock). The numerical model combines uncoupled, linear pressure diffusion with the Mohr‐Coulomb failure criterion to simulate IIS. We demonstrate that connected pathways of relatively high hydraulic diffusivity in heterogeneous media (particularly in fractured rock domains) allow the threshold triggering pressure to propagate more rapidly than predicted by the effective hydraulic diffusivity. As a result, the seismic diffusivity is greater than the effective hydraulic diffusivity in heterogeneous porous media, possibly by an order of magnitude or more. Additionally, we present a case study of IIS near Soultz‐sous‐Forêts where seismic diffusivity is found to be at least one order of magnitude larger than the effective hydraulic diffusivity.

It is well established that the propagation of seismicity and the triggering front (as observed in r-t plots) is approximately described by a diffusion-like process (Shapiro, 2015;Shapiro et al., 1997;Talwani & Acree, 1985). As such, expressions describing the diffusive propagation of the triggering front can be fit to the upper envelope of the seismic data cluster, producing an estimate of diffusivity (Delepine et al., 2004;Hummel & Müller, 2009;Hummel & Shapiro, 2012Rothert & Shapiro, 2003;Segall & Lu, 2015;Shapiro, 2015;Shapiro & Dinske, 2009a, 2009bShapiro & Müller, 1999;Shapiro et al., 1997Shapiro et al., , 2002. For clarity, we refer to this estimate of diffusivity as the "seismic diffusivity," borrowing the term from Talwani and Acree (1985). There, the authors suggested that estimates of seismic diffusivity (associated with the spatiotemporal patterns of IIS) are in fact accurate estimates of the effective hydraulic diffusivity of the injection formation. If true, this indicates that methods for estimating seismic diffusivity are in fact seismicity-based approaches for subsurface hydraulic characterization. Since Talwani and Acree (1985), numerous other studies have also implicitly assumed that seismic diffusivity and effective hydraulic diffusivity are equivalent, and thus the term "seismic diffusivity" is not widely used Haagenson et al., 2018).
An important question that arises is whether the aforementioned equivalence between seismic diffusivity and effective hydraulic diffusivity holds in a heterogeneous porous medium. When heterogeneity is represented as a spatially correlated random field with a well-defined correlation length, it is generally accepted that at scales much larger than the correlation length, the behavior of the mean pressure and fluid flux fields can be described by effective permeability and hydraulic diffusivity tensors. Similar frameworks exist for defining block-scale effective conductivities for large grid blocks in field-scale flow models based on the underlying heterogeneity structure or fracture network topology (Botros et al., 2008;Long et al., 1982;Renard & De Marsily, 1997;Sanchez-Vila et al., 2006;Wen & Gómez-Hernández, 1996;Zimmerman & Bodvarsson, 1996a). Most field-scale modeling of fluid flow in hydrogeology implicitly employs equivalent homogeneous representations within geological units and only represents material property variations across distinct geological units explicitly. Moreover, field-scale hydraulic tests are commonly interpreted based on analytical solutions for homogeneous media and the properties estimated from these tests are assumed to represent effective properties at the test scale (Gelhar, 1993). In radial flow, theoretical and computational analyses in heterogeneous media (e.g., Guadagnini et al., 2003;Naff, 1991) demonstrate that apparent conductivities approach constant values equal to the effective conductivities within a few correlation lengths from the well. Analytical and numerical modeling studies of IIS and frameworks for interpreting the spatiotemporal patterns of IIS implicitly adopt equivalent homogeneous representations of hydraulic properties by relying on the assumed equivalence between the seismic diffusivity and effective hydraulic diffusivity of the subsurface (Brown et al., 2017;Catalli et al., 2016;Dempsey & Riffault, 2019;Keranen et al., 2014;Langenbruch & Zoback, 2016;Langenbruch et al., 2018;Pollyea et al., 2019;Riffault et al., 2018;Segall & Lu, 2015;Shapiro, 2015;Shapiro & Dinske, 2009a, 2009bShapiro & Müller, 1999;Shapiro et al., 1997).
Our specific goal in this paper is to demonstrate that in highly heterogeneous porous media, there is a clear distinction between the seismic diffusivity and the effective hydraulic diffusivity. In reality, the triggering front is not radially symmetric in heterogeneous domains. This is because pressure increments are expected to propagate preferentially through pathways of relatively high hydraulic diffusivity in heterogeneous media. As a result, the farthest radial distance to which the threshold triggering pressure increment propagates (which is synonymous with the triggering front in the context of r-t plots) will be farther than the radial distance at which the directionally averaged pressure increment equals p t (which propagates according to the effective hydraulic diffusivity). Therefore, the triggering front observed in r-t plots will propagate more rapidly than predicted by the effective hydraulic diffusivity. Put another way, seismicity will propagate more rapidly in a realistic, heterogeneous domain than it would in an equivalent homogeneous domain that employs the effective hydraulic diffusivity, since the latter would inevitably underestimate the farthest radial distance reached by the threshold triggering pressure increment (p t ) at any given time.
Using numerical simulations of fluid flow and induced seismicity, we investigate this possible distinction, the degree to which these two quantities may differ, and what subsurface conditions may exacerbate the difference. After providing a description of the computational framework in Section 2, we present results of simulated fluid injection and induced seismicity in two different types of heterogeneous domains: smoothly varying fields of random hydraulic diffusivity derived using Sequential Gaussian Simulation (SGS) (i.e., SGS domains described in Section 2.1) and three-dimensional discrete fracture networks in rock matrix (i.e., discrete fracture network and matrix [DFNM] domains described in Section 2.2). In Section 4, we present a case study of the fluid injection experiment near Soultz-sous-Forêts, France from the year 2000, which further illustrates the potential distinction between the seismic and effective hydraulic diffusivity (Dezayes et al., 2010;Dorbath et al., 2009;Genter et al., 2010;Meller & Ledésert, 2017).

Computational Framework
The correlated random fields of hydraulic diffusivity are generated using a SGS algorithm, producing smoothly varying fields (which are a common but relatively simple representation of heterogeneity) (Müller & Schüler, 2020). Domains of highly fractured rock (generated using dfnWorks ) are more realistic representations of the kind of formation in which IIS typically occurs. We refer to these HAAGENSON AND RAJARAM 10.1029/2021JB021768 3 of 25 domains as a DFNM domains, because we explicitly model fluid flow in both the fracture network and surrounding rock matrix (Birdsell et al., 2018;Haagenson et al., 2018).
For simplicity, we consider only heterogeneous media with large-scale, isotropic effective permeability and diffusivity tensors (i.e., statistically isotropic correlation structures in the SGS domains and uniform orientation distributions for fractures in the DFNM domains). These models are readily generalizable to represent large-scale anisotropy resulting from either anisotropic spatial correlation in the SGS domains or non-uniform orientation distributions (e.g., families of similarly aligned fractures) in the DFNM domains. However, our purpose here is mainly to illustrate the influence of heterogeneity on the propagation of the triggering front, which is expected to occur with or without anisotropy.
We model fluid flow through a porous medium in response to a point source of fluid injection based on the classical pressure diffusion equation (Bear, 1972;Charbeneau, 2006;Cleary, 1977;De Marsily, 1986;Rice & Cleary, 1976;Rudnicki, 1986;Verruijt, 2013), given as where ρ f is the fluid density, ϕ is the porosity of the porous medium, β f and β m are the compressibilities of the fluid and porous medium respectively, p is the increment of pore fluid pressure (with respect to the initial, static pressure distribution), κ is the intrinsic permeability of the porous medium, μ is the dynamic viscosity of the fluid, and Q m is the source or sink of fluid mass. Although some of these parameters (particularly permeability) can exhibit pressure dependence, we use constant parameter values such that Equation 1 becomes a linear diffusion equation. Our goal is to highlight the influence of heterogeneity on patterns of IIS within the framework of linear diffusion. We intend to investigate the nuanced behavior of nonlinear fluid flow and IIS in future work. In this form, Equation 1 represents a diffusion equation of pressure increment, where the well-known hydraulic diffusivity (Bear, 1972;De Marsily, 1986) is defined as where b is the fracture aperture, Q f is the source or sink of fluid mass per unit area of the fracture, and L m is the fluid leak-off rate per unit area from the fracture into the surrounding rock matrix (Chaudhuri et al., 2013;Murphy et al., 2004). Equation 3 employs the well-known local cubic law for fracture transmissivity (Zimmerman & Bodvarsson, 1996b). Although there are limitations to the local cubic law at smaller scales, it is widely used in discrete fracture network (DFN) models (Adler et al., 2013;Frampton et al., 2019; with a constant aperture (b) within individual fractures (which may also be interpreted as the equivalent hydraulic aperture). In Section 2.2.2 below, we further describe how Equation 3 can be expressed in the form of a diffusion equation and solved in combination with the pressure diffusion equation in the rock matrix given in Equation 1.
For simulations of fluid injection, we have developed a numerical model to solve Equations 1 and 3 using FEniCS-a general-purpose, open-source finite element method software (Alnaes et al., 2015), which has been previously applied to a wide range of geoscientific problems including subsurface fluid flow and generalized poroelasticity . Each domain used in this study is cubic, with sides measuring 2 km in length. Fluid is injected at the center vertex of the numerical mesh (as a point source) at a constant rate (Q) of 25 m 3 per hr for a period (t e ) of 12 hrs. These parameters were selected to reflect a realistic scenario of fluid injection (Shapiro & Dinske, 2009b;Shapiro et al., 1997Shapiro et al., , 2005, while also minimizing potential boundary effects by employing a sufficiently large domain (with domain boundaries located a distance of at least 12 eff e D t away from the injection location, where D eff is the effective hydraulic diffusivity of the heterogeneous domain). The initial condition and all boundary conditions are set to a hydrostatic pressure field. A summary of the parameters used in the numerical simulations are given in Table 1.
Following the approach of Rinaldi and Nespoli (2017), Catalli et al. (2016), and Shapiro (2015), we track seismicity using a set of weak points, which are seeded randomly within 1 km of the injection location. Each weak point represents a potential location for seismicity. In the SGS domain, each weak point is randomly assigned a strike and dip angle (to represent a hypothetical fracture at that location), whereas the weak points in the DFNM domain are seeded exclusively along fractures and are assigned the corresponding fracture's strike and dip angle. The weak point triggering pressure increment ( wp t p ), which is the increase in fluid pressure above the initial fluid pressure that would trigger a seismic event, is found using the well-known Mohr-Coulomb failure criterion: where f is the coefficient of static friction of the fracture, p i is the initial fluid pressure at the weak point, and σ n and τ are the normal compressive stress and shear stress acting on the fracture respectively. Note that the term    ( ) Equation 4 is the effective stress acting normal to the fracture plane at the weak point location. A decrease in the effective normal stress due to an increase in the fluid pressure is commonly assumed as the dominant physical mechanism that triggers seismicity (Dempsey & Riffault, 2019;Healy et al., 1968;Keranen et al., 2014;Riffault et al., 2018). The coefficient of friction is allowed to vary between 0.6 and 0.7 to capture a realistic range of values (Jaeger, 1959;Talwani & Acree, 1985). The initial pressure profile is assumed to be hydrostatic, and is based on a constant fluid density (ρ f = 998 kg/m 3 ) and a water table located 100 m below the ground surface. The depth of the fluid injection (H inj ) is stipulated as 4 km. We have assumed that the local stress tensor is defined by a normal faulting, lithostatic condition where the vertical stress component is defined as σ v = ρ s g(H inj − z) and both horizontal components as σ h = 0.6σ v (Zoback, 2010). The rock density of the overburden (ρ s ) is assumed to be 2,300 kg/m 3 . More details on the seeding of weak points in the SGS and DFNM domains are provided in Text S1.
Pressure increments associated with induced seismicity are commonly assumed to lie in the range of 0.01 and 0.1 MPa. A pressure increment less than 0.01 MPa is either not significant enough to trigger seismicity or is often met by processes other than fluid injection (such as background fluid flow or tidal forcing), whereas a pressure increment larger than 0.1 MPa is likely to occur on fractures or faults with orientations that are not favorable for shear slip (Dempsey & Riffault, 2019;Gischig & Wiemer, 2013;Goebel et al., 2017;Keranen et al., 2014;Shapiro, 2015). Thus, weak points with a triggering pressure increment outside this range are removed from the final set of weak points used in the simulation. It is clear then, that the threshold triggering pressure increment associated with the triggering front will in fact be the minimum of the range of weak point triggering pressure increment values (i.e.,   ( ) 0.01 MPa wp t t p min p ).
To track seismicity, each weak point is evaluated for potential failure in each time step of the simulation. This approach results in a synthetic IIS data set, which allows us to produce r-t plots in order to investigate the spatiotemporal patterns of seismicity and estimate the seismic diffusivity of the heterogeneous domain (following the approach described in Section 2.4).

SGS Domains
There is a long history of previous work on flow and transport in heterogeneous porous media based on numerical simulations in computer-generated spatially correlated random fields (Deutsch & Journel, 1998;Tompson et al., 1989). To generate such a field for hydraulic diffusivity, we first generate a logarithm of  permeability field (ln κ, where κ has units of m 2 ) with an isotropic, Gaussian covariance function using the GeoStat-Framework (i.e., GSTools module in Python) (Müller & Schüler, 2020). We simulate 10 realizations of the SGS domain at four different values of the variance of ln κ (Var(ln κ) = 0.53, 1.325, 2.65, and 5.3), for a total of 40 SGS domain realizations. An example of the resulting ln κ field is shown in Figure 2. A correlated, random field of permeability (κ) for use in Equations 1 or 2 is readily obtained by taking the exponential of the ln κ field.
After generating the field of permeability, the correlated, random field of hydraulic diffusivity is found using the well-established definition of hydraulic diffusivity given in Equation 2. When evaluating this expression, we have assumed typical, constant values for μ (8.9 × 10 −4 Pa⋅s), ϕ (0.05), β f (4.4 × 10 −10 Pa −1 ), and β m (7.0 × 10 −11 Pa −1 ) that are reflective of a highly fractured, low permeability rock in which IIS most often occurs (Bear, 1972;Charbeneau, 2006;De Marsily, 1986;Freeze & Cherry, 1987). Note that permeability can vary spatially by several orders of magnitude within a heterogeneous domain, whereas porosity varies within a relatively narrow range in a given formation (Gelhar, 1993). Therefore, using a spatially variable permeability along with a uniform porosity is an acceptable approximation that is common in studies of fluid flow in heterogeneous porous media. A summary of the model parameters specific to the SGS domains are given in Table 2.

DFNM Domains
In a domain of highly fractured, low permeability rock, the fracture network represents a well-connected pathway of relatively high hydraulic diffusivity through which pressure increments can propagate rapidly. Since most cases of IIS occur in fractured rock, it is critical to consider models of subsurface heterogeneity that explicitly represent these connected fracture networks. In the following sections, we describe the methodology specifically employed for the fractured rock domain.

Fracture Network and Numerical Mesh Generation
The DFNM domain is fully three-dimensional, with the fracture network represented by penny-shaped fractures of finite thickness embedded within the rock matrix. This methodology has been previously presented (Birdsell et al., 2018;Haagenson et al., 2018) and is also referred to as an upscaled discrete fracture-matrix model by Sweeney et al. (2020). The primary benefits of this approach, as opposed to modeling the fracture network alone, is that we can directly account for fluid flow in the rock matrix (including leak-off from the fracture network) while keeping the mathematical and numerical formulations relatively simple. Interface finite element formulations can handle narrow fractures treated as discontinuities (Abushaikha et al., 2015;Berre et al., 2019;Fumagalli & Scotti, 2013;Geiger et al., 2010;Odsaeter et al., 2019); however, these approaches typically limit the domain to either two dimensions or simple DFN geometries in three dimensions (Sweeney et al., 2020). The DFNM domain used in this study is HAAGENSON AND RAJARAM 10.1029/2021JB021768 6 of 25   Figure 3. More details on the generation of the DFNM domain are provided in Text S2. A summary of the model parameters specific to the DFNM domain is presented in Table 3.
The zones of relatively high hydraulic diffusivity in the SGS realizations with Var(ln κ) = 5.3 (i.e., the largest variance considered) have hydraulic diffusivity values comparable to those found in the fractures of the DFNM domain. The major difference between the domains is that zones of high hydraulic diffusivity are well connected in the DFNM domain.
To illustrate this, we plot the distribution of hydraulic diffusivity within the fracture network of the DFNM domain in Figure 4b, and compare it with the locations of equivalently high hydraulic diffusivity in an example of the SGS domains (with Var(ln κ) = 5.3) in Figure 4a. Notice that zones of relatively high hydraulic diffusivity are not well connected in the SGS case, whereas they are in the DFNM domain. We expect that the enhanced connectivity in the DFNM domain will cause pressure increments to diffuse rapidly within the fracture network, leading to more rapid propagation of seismicity compared to the SGS domains.

Flow and Seismicity in a Fracture Network
Fluid flow in the rock matrix portion of the DFNM domain is governed by Equation 1. Here, we discuss the treatment of Equation 3 for modeling fluid flow in the fracture network. The DFNM domain contains both the fracture network and rock matrix (as discussed in Section 2.2.1), allowing us to drop L m from Equation 3 as leak-off from the fracture network is automatically considered by allowing fluxes across fracture-matrix interfaces. The storage term in Equation 3 may be expanded following the approach of Rutqvist et al. (1998), giving where K n is the normal stiffness of the fracture, which can be measured using laboratory experiments (Bandis et al., 1983). Following the approach of Chaudhuri et al. (2013) and others (Birdsell et al., 2015(Birdsell et al., , 2018Bower & Zyvoloski, 1997;Pandey & Rajaram, 2016;Pandey et al., 2017), we define the permeability of the fractures as    where b p is the model fracture cell width. The numerical mesh of the DFNM domain is recursively refined around the location of the fractures, giving the model fractures a final thickness of b p . This process is described in more detail in Text S2. Typically, fracture permeability is given by the well-known expression κ frac = b 2 /12; however, scaling this expression with b/b p in Equation 6 accurately accounts for the discrepancy between b and b p . This approach in fact rigorously ensures that the model fractures are hydraulically equivalent to the actual fractures (Birdsell et al., 2015(Birdsell et al., , 2018Bower & Zyvoloski, 1997;Chaudhuri et al., 2013;Pandey & Rajaram, 2016;Pandey et al., 2017;Sweeney et al., 2020). Substituting Equation 6 into Equation 5 and rearranging gives This is the final form of the governing equation for pressure diffusion in the fracture cells of the DFNM domain, where Q v is the source or sink term of fluid volume per unit volume of the fracture. Similar to Equation 1, we assume that the parameters in Equation 7 are independent of fluid pressure leading to a linear diffusion equation. Some fracture parameters (particularly fracture aperture) are known to be pressure and stress dependent (e.g., normal dilation [Bandis et al., 1983] or shear dilation [Rong et al., 2016;Ye & Ghassemi, 2018]), which may influence the behavior of fluid flow and propagation of seismicity in the DFNM domain. Since the main objective here is to investigate the influence of heterogeneity, we use a linear form of Equation 7. We intend to investigate the nuanced influence of nonlinear fluid flow in fracture networks on the behavior of seismicity in future work. From Equations 6 and 7, we see that the hydraulic diffusivity of the fractures can be readily defined as

Estimating Effective Hydraulic Properties of Heterogeneous Domains
The effective permeability and hydraulic diffusivity of the heterogeneous domains were determined using numerical permeameter and injection tests, which are described in detail in Text S3 (Brezzi & Fortin, 2012;Brezzi et al., 1985;Chrysikopoulos, 1995;Dagan, 1979Dagan, , 1989Dykaar & Kitanidis, 1992a, 1992bFernàndez-Garcia et al., 2005;Gelhar, 1993;Gutjahr et al., 1978;Paleologos et al., 1996;Traverso et al., 2013;Younes et al., 2010;Zhang, 2001). In short, the effective permeability of each heterogeneous domain (κ eff ) was rigorously estimated using a numerical permeameter test in each axial direction and then averaged to find a single value for the entire domain (i.e., κ eff = (κ xx + κ yy + κ zz )/3, based on the expectation of isotropic HAAGENSON AND RAJARAM 10.1029/2021JB021768 8 of 25 macro-scale hydraulic properties). The effective hydraulic diffusivity (D eff ) of each SGS domain was then calculated using Equation 2, with the effective permeability value (κ eff ) from the numerical permeameter test and the uniformly applied values of all other parameters. We confirmed these estimated values of effective permeability and hydraulic diffusivity using a numerical injection test, comparing the directionally average fluid pressure profile with an analytical solution for the fluid pressure profile in response to a point source of fluid mass in a isotropic and homogeneous medium. This analytical solution is given by where erfc is the complementary error function (Carslaw & Jaeger, 1959). Since the DFNM domain contains widely disparate values of hydraulic and mechanical properties between the rock matrix and fractures, we cannot evaluate Equation 2 directly as we have done for the SGS domains. Therefore, the effective hydraulic diffusivity of the DFNM domain is estimated by fitting Equation 9 to the profile of directionally averaged fluid pressure during a fluid injection simulation using a nonlinear regression algorithm, as described in Text S3.

Estimating Seismic Diffusivity
The seismic diffusivity is typically estimated using expressions that track the location of the triggering front during fluid injection. For example, the expression   4 t r Dt (where D is the diffusivity value being estimated) is very commonly used (Antonioli et al., 2005;Chen et al., 2012;Goebel et al., 2017;Goertz-Allmann et al., 2017;Haffener et al., 2018;Hummel & Shapiro, 2012;Ingebritsen & Manning, 2010;Shapiro & Dinske, 2009a;Yong et al., 2018;Yu et al., 2019). However, it is more appropriate to use the analytical solution to linear pressure diffusion in a three dimensions in response to a point source of fluid since this approach incorporates physically meaningful parameters such as the fluid injection rate. Assuming that the triggering front is equivalent to the isobaric surface with pressure equal to the threshold triggering pressure increment (p t ), then the arrival time of the triggering front at any given radius is given by where D s is the seismic diffusivity associated with the diffusive propagation of the triggering front and erfc −1 is the inverse of the complementary error function (Carslaw & Jaeger, 1959). Equation 10 may be fit to the upper envelope of the seismic data cluster found in r-t plots in order to produce an estimate of the seismic diffusivity. It is clear from Equation 10 that the seismic diffusivity can in fact be interpreted as the hydraulic diffusivity of a hypothetical, homogeneous domain that (when subjected to the same injection conditions) approximately produces the same triggering front behavior as observed in the r-t plot based on seismic data from the real-world, heterogeneous domain. This interpretation agrees well with the concept of seismic diffusivity presented by Talwani and Acree (1985) and with the methods of Shapiro et al. (1997) and Shapiro (2015). More details on estimating the seismic diffusivity are provided in Text S4.

Results of Pressure Diffusion and Induced Seismicity
The results from the numerical simulations of both the SGS and DFNM domains are presented in the following two sections. We present results of both pressure diffusion and induced seismicity in three-dimensional domains, which illustrates the general influence of heterogeneity. We also examine the behavior of the simulated IIS in r-t plots, which provides insights into the concept of seismic diffusivity and its potential distinction from the effective hydraulic diffusivity of the heterogeneous domains.

SGS Domain Results
We begin by investigating the influence of heterogeneity on the patterns of IIS in the context of SGS domains, to gain a fundamental understanding of the influence of heterogeneity on the spatiotemporal patterns of IIS in relatively simple representations of subsurface heterogeneity. As already mentioned, we performed fluid injection simulations in 40 different realizations of the SGS domain, with 10 realizations at each value of Var(ln κ) considered (i.e., 0.53, 1.325, 2.65, 5.3). Figure 5 shows the pressure solution from one example realization of the SGS domain at each value of Var(ln κ) to illustrate the three-dimensional patterns of pressure diffusion. Cutaway sections show the internal distribution of fluid pressure. Portions of the domain with pressure increment below the threshold triggering pressure increment (p t ) are not shown. In essence, the outer, dark blue contour (which represents a pressure increment value of p t ) shows the location of the triggering front in three-dimensional space, revealing the clear influence of heterogeneity. Figure 5 shows that as the value of Var(ln κ) increases, the pressure contours and thus the triggering front deviate farther from radial symmetry. We see in the Var(ln κ) = 5.3 case shown in Figure 5d that p t has propagated much farther in certain directions than others, which correspond to zones of relatively high hydraulic diffusivity. Figure 6 shows the r-t plot from one of the SGS domain realizations with Var(ln κ) = 5.3, overlaid with the maximum radial distance reached by the threshold triggering pressure increment (denoted by r max (p = p t )) and the radial distance at which the directionally averaged pressure increment (i.e., the averaged pressure increment over a sphere of a given radius) equals the threshold triggering pressure increment (denoted by r(p avg = p t )). First, we see clearly that the farthest occurring seismicity at any given time approximately coincides with r max (p = p t ). This is the expected result, since weak points with triggering pressure increments HAAGENSON AND RAJARAM 10.1029/2021JB021768 10 of 25 near the threshold value of p t will trigger soon after the arrival of p t . Comparing the locations of r max (p = p t ) and r(p avg = p t ), we clearly see that r max (p = p t ) propagates much more rapidly. This is due to the presence of zones of relatively high hydraulic diffusivity that allow pressure increments to diffuse rapidly in certain directions away from the injection location. This rapid propagation of p t will directly influence r max (p = p t ) while the value of r(p avg = p t ) lags behind (due to its very nature as a directional average, meaning it is influenced by the entire distribution of hydraulic diffusivity and not just the relatively high values). . Example r-t plot from an Sequential Gaussian Simulation realization with Var(ln κ) = 5.3. The plot is overlaid with the maximum radial distance to which the threshold triggering pressure increment has propagated at any given time (i.e., r max (p = p t ) in red) and the radial distance at which the directionally averaged pressure increment equals the threshold triggering pressure increment (i.e., r(p avg = p t ) in blue).

Figure 7.
Left: Plots of the radial distance at which the directionally averaged pressure increment equals the threshold triggering pressure increment (i.e., r(p avg = p t )) and the maximum radial distance reached by the triggering pressure at any given time (i.e., r max (p = p t )). Light blue and red lines show results from each of 10 Sequential Gaussian Simulation domain realizations. Dark blue and red lines show ensemble mean. Right: Plots of the ensemble means of r(p avg = p t ) and r max (p = p t ). Two envelopes about the ensemble mean are plotted: one standard deviation (dark shading) and two standard deviations (dark shading).
r(p avg = p t ) across the 10 realizations; although, there is more variability in the behavior of r max (p = p t ). This is expected, as r max (p = p t ) estimates a maximum radial distance of a certain pressure increment value, whereas r(p avg = p t ) is inherently a directionally averaged variable and is, therefore, less sensitive to rapid propagation of p t in a particular direction. If p t propagates rapidly in a certain direction relative to others, that will be directly reflected in r max (p = p t ), but will not heavily influence r(p avg = p t ).
Despite the increase in variability across realizations, it is clear that the deviation of r max (p = p t ) from r(p avg = p t ) increases systematically as the degree of heterogeneity increases (i.e., as Var(ln κ) increases).
From the 10 realizations at each value of Var(ln κ), we also calculated the standard deviations of r max (p = p t ) and r(p avg = p t ), representing their variability across realizations. In Figure 7 (right), we show the ensemble mean of both r max (p = p t ) and r(p avg = p t ) along with two envelopes: one standard deviation above or below the ensemble mean and two standard deviations above or below the ensemble mean. The envelopes around r(p avg = p t ) are much narrower than those around r max (p = p t ), indicating that the propagation of r max (p = p t ) may be more sensitive to the influence of heterogeneity. For the lower values of Var(ln κ) (i.e., 0.53 and 1.325 where the domains are relatively homogeneous), the envelopes show a significant amount of overlap, indicating that they are not statistically distinct from one another. However, for higher values of Var(ln κ) (i.e., 2.65 and 5.3 where the domains are relatively heterogeneous), the envelopes do not overlap at late time during the simulations, indicating that r max (p = p t ) and r(p avg = p t ) are statistically distinct from each other, even in the sense of ensemble means.
As commonly done in IIS literature, we use the spatiotemporal patterns of IIS and the propagation of the triggering front to estimate the seismic diffusivity of the heterogeneous domain. The details of our approach are outlined in Section 2.4. To exemplify this approach, we fit Equation 10 to the upper envelope of the seismic data cluster (as shown by the dashed red line) in Figure 8a using the same example r-t plot shown in Figure 6. We see that the triggering front (associated with the seismic diffusivity) generally reflects the approximate behavior of the maximal radial distance to which the threshold triggering pressure increment has propagated (i.e., r t ≈ r max (p = p t ), where r t is the location of the triggering front according to Equation 10).
Next, we consider the propagation of the triggering front in a hypothetical homogeneous domain with uniform diffusivity equal to the effective hydraulic diffusivity of the SGS domain, which is overlaid in Figure 8a as a dashed blue line. We clearly see that this would drastically underestimate the actual location of the triggering front, the maximal radial distance to which p t propagates, and the farthest reaches of seismicity at any given time. In fact, Equation 10 with the effective hydraulic diffusivity closely coincides with r(p avg = p t ). This is consistent with the concept of the effective hydraulic diffusivity, as it is associated with the evolution of the propagation of directionally averaged pressure contours of any pressure increment value (including the contour of p t ).
We estimate the seismic diffusivity of each SGS domain, and plot these values against Var(ln κ) in Figure 8b. Clearly, the seismic diffusivity of the domain increases with Var(ln κ). This indicates that the rate at which seismicity propagates will also increase with Var(ln κ). However, we expect the effective hydraulic diffusivity of the SGS domains to increase with Var(ln κ) as well. We therefore evaluate the ratio of D s /D eff for the SGS domains, plotted in Figure 8c. It is clear that, while the variability of this ratio across the SGS realizations increases with Var(ln κ), the ensemble average of this ratio systematically increases with Var(ln κ) as well.
HAAGENSON AND RAJARAM 10.1029/2021JB021768 12 of 25 For any value of Var(ln κ), there are some realizations in which the regions of relatively high hydraulic diffusivity are not well connected. In such realizations, the ratio D s /D eff will be only slightly above one, as pathways of high hydraulic diffusivity through which pressure increments can rapidly propagate are limited. However, in realizations where regions of relatively high hydraulic diffusivity are well connected, the ratio D s /D eff can be quite large and will increase with Var(ln κ). Overall, the ensemble average indicates a systematic increase of the ratio D s /D eff with Var(ln κ).

DFNM Domain Results
Next, we consider a DFNM domain, which is a better representation of fractured rock formations in which IIS commonly occurs. In the study of SGS domains, we mainly focused on showing that the seismic diffusivity (which is associated with the propagation of r max (p = p t )) and the effective hydraulic diffusivity (which is associated with the propagation of r(p avg = p t )) of the highly heterogeneous SGS domains are not equivalent. In the DFNM domain, the zones of relatively high hydraulic diffusivity (i.e., the fractures) are well connected (see Figure 4). This develops well-connected pathways of high hydraulic diffusivity (which were only occasionally present in the SGS domains) through which the pressure increment from the injection location can rapidly propagate, causing seismicity to spread rapidly as well. Hence, the discrepancy between the seismic diffusivity and the effective hydraulic diffusivity of the DFNM domain may be larger than in the SGS cases. To evaluate this, we calculate the ratio D s /D eff for the DFNM domain, and compare it to the values found in the study of the SGS domains.
Additionally, we highlight one major implication of the discrepancy between the seismic diffusivity and the effective hydraulic diffusivity. If the spatiotemporal patterns of seismicity are not in fact accurately described by the effective hydraulic diffusivity of a realistic, heterogeneous domain, then the common approach of using homogeneous domains (or layers of homogeneous domains) in numerical simulations of induced seismicity (Brown et al., 2017;Catalli et al., 2016;Dempsey & Riffault, 2019;Keranen et al., 2014;Langenbruch & Zoback, 2016;Langenbruch et al., 2018;Pollyea et al., 2019;Riffault et al., 2018) may be limited. To investigate this, we compare the results of simulated fluid injection in two domains: the DFNM domain and an equivalent homogeneous domain with uniform hydraulic diffusivity equal to the effective hydraulic diffusivity of the DFNM domain. As mentioned in Section 2.3, the effective hydraulic diffusivity of the DFNM domain was estimated using numerical permeameter and fluid injection tests and was found to exhibit isotropic behavior. Thus, the equivalent homogeneous domain used here is isotropic as well.
The results of pressure diffusion in response to fluid injection for both domains are shown in Figure 9. It is clear that the threshold triggering pressure increment has propagated to a significantly greater distance away from the injection location in the DFNM case. This is the expected result since fluid flows radially outward in the homogeneous domain, whereas fluid flows primarily in the fracture network of the DFNM domain. Forcing the stipulated injection rate into the relatively small volume of the of the fracture network requires that pressure increments from the injection location propagate more rapidly in the DFNM case than in the homogeneous one. Hence, the threshold triggering pressure increment has propagated to a much farther distance by the end of the injection simulation in the DFNM case.
In Figure 10, we show the locations of induced seismicity triggered during the numerical simulation of fluid injection for each domain. Since the threshold triggering pressure increment propagates much farther in the DFNM case, seismicity in the DFNM case occurs at greater distances away from the injection location as well. Additionally, we see two unique behaviors in the DFNM case. First, events associated with a relatively high triggering pressure increment were triggered much more often in the DFNM case. Since the stipulated fluid injection rate must be accommodated almost entirely by the DFN (which has a relatively small volume), large pressure increments from the injection location must propagate farther into the DFNM domain than in the homogeneous one. This triggers weak points in the simulation of the DFNM domain that are less prone to fail (due to a higher coefficient of friction or an unfavorable fracture orientation that hinders slip). Second, we see in bottom left corner of Figure 10e that some fractures of the DFNM domain have actually been pressurized to the threshold triggering pressure increment (p t ) while not experiencing any seismicity. The unfavorable orientation of these fractures limits the potential for seismicity. Still, these fractures act as efficient pathways to facilitate pressure increment propagation. Fluid flow through these fractures may therefore lead to seismicity along a connected fracture that is more favorably oriented for  . The spatial extent to which threshold triggering pressure increment (p t ) extends (which corresponds to the plotted portions in Figure 9) is shaded in dark gray. failure, as we see in Figure 10f. These are examples of the more subtle influences of fracture networks on the behavior of IIS.
If we look at the synthetic data set of seismicity from both simulations as an r-t plot (shown in Figure 11), we clearly see very different behavior between the equivalent homogeneous and DFNM domains. As expected, seismicity has propagated much farther in the DFNM domain than in the equivalent homogeneous domain. This suggests that equivalent homogeneous domains (which are commonly used in studies of IIS [Brown et al., 2017;Catalli et al., 2016;Dempsey & Riffault, 2019;Keranen et al., 2014;Langenbruch & Zoback, 2016;Langenbruch et al., 2018;Pollyea et al., 2019;Riffault et al., 2018]) may be limited in their ability to replicate the spatiotemporal patterns of IIS as observed in r-t plots from real-world or heterogeneous domains. Next, we examine the behavior of the triggering front in each domain, using Equation 10 as described in Section 2.4. As expected, the triggering front in the equivalent homogeneous domain is accurately described using the effective hydraulic diffusivity. In this case, the triggering front is radially symmetric and thus there should not be any distinction between the seismic diffusivity and the effective hydraulic diffusivity. In the DFNM domain, however, we see that the location of the triggering front is vastly underestimated by Equation 10 using the effective hydraulic diffusivity of the domain, as was the case in the highly heterogeneous SGS domains. This again suggests that the seismic diffusivity associated with the diffusive propagation of the triggering front (which was determined to be 3.0 m 2 /s in the DFNM domain) is larger than the effective hydraulic diffusivity (which was determined to be 0.29 m 2 /s in the DFNM domain).
Finally, we consider the ratio of D s /D eff as we did in the study of SGS domains (see Figure 8c for the comparison of this ratio from the SGS domains), which produces a value of 10.2 in the DFNM domain. This shows that in the case of fractured rock, the distinction between the seismic and effective hydraulic diffusivity can be up to one order of magnitude or more. This is considerably larger than the largest ratio found in the SGS cases (i.e., D s /D eff = 3.4), which indicates that the well-connected pathways of the fracture network allow pressure increments to propagate more rapidly in the DFNM domain than in any of the SGS domains considered in this study. The ratio of D s /D eff for the DFNM domain is just an illustrative example and higher ratios may be possible for other fracture network topologies or densities. A simple sensitivity analysis to gain insight into the potential uncertainty associated with the estimate of seismic diffusivity of the DFNM domain is presented in Text S5.

Case Study: Soultz-Sous-Forêts
Numerous studies of IIS have previously indicated that tracking the location of the triggering front as observed in r-t plots may produce reasonable estimates of the effective hydraulic diffusivity at the injection site, mostly using the approach described by Shapiro et al. (1997) (Antonioli et al., 2005;Chen et al., 2012;Delepine et al., 2004;Goebel et al., 2017;Goertz-Allmann et al., 2017;Haffener et al., 2018;Hummel & Müller, 2009;Hummel & Shapiro, 2012Improta et al., 2015;Ingebritsen & Manning, 2010;Rothert & Shapiro, 2003;Segall & Lu, 2015;Shapiro & Dinske, 2009a, 2009bShapiro & Müller, 1999;Shapiro et al., 2002;Yong et al., 2018;Yu et al., 2019). Although these previous studies have implicitly assumed that the seismic diffusivity is essentially equivalent to the effective hydraulic diffusivity, the numerical study presented here has clearly shown that these two quantities are likely distinct in heterogeneous domains. To further illustrate this potential distinction, we analyze the seismic data set collected during the fluid injection experiment near Soultz-sous-Forêts, France in the year 2000.
The lower reservoir of the enhanced geothermal system (EGS) site near Soultz-sous-Forêts, France consists of three wells (called GPK2, GPK3, and GPK4) located in a highly fractured, crystalline rock formation HAAGENSON AND RAJARAM 10.1029/2021JB021768 16 of 25 approximately 4,000-5,000 m below the ground surface (Dezayes et al., 2010;Genter et al., 2010;Meller & Ledésert, 2017). The site was selected for EGS due to the presence of a thermal anomaly, with downhole temperatures reaching 200 °C in each of the three wells (Meller & Ledésert, 2017). Hydraulic stimulation was performed at each well at different times during the construction of the EGS site: GPK2 in 2000, GPK3 in 2003, and GPK4 in 2004(Dorbath et al., 2009Meller & Ledésert, 2017). A full description of the site geology and wellbore details are given by Meller and Ledésert (2017) and Dezayes et al. (2010). Fluid injection and seismicity data have been presented by Genter et al. (2010) and Cuenot et al. (2008), in addition to being freely available through the EPOS-IP Anthropogenic Hazards data repository (Leptokaropoulos et al., 2019). Here, we present a case study of the fluid injection experiment performed at the GPK2 well in the year 2000.
Water was injected at an average rate of 44.6 L per second (see Figure 12a) for a period of approximately 6 days. The open hole interval of the injection well ranged from 4,400 to 5,100 m below the ground surface. Using numerous surface and downhole seismometers, seismicity was tracked in space and time during the fluid injection experiment, which reached up to 1 km away from the injection location. During just the first day of fluid injection, seismicity occurred as far as 700 m away from the injection location. The r-t plot of this seismic data set was shown earlier in Figure 1.
In order to characterize the hydraulic diffusivity of the EGS reservoir, Delepine et al. (2004) previously studied the spatiotemporal patterns of IIS found in this seismic data set using the expression which was originally suggested by Shapiro et al. (1997). By fitting Equation 11 to the upper envelope of the seismic data cluster shown in Figure 12, the authors estimated the seismic diffusivity to be 0.15 m 2 /s. (Note that the authors did not use the term "seismic diffusivity," but instead assume that the spatiotemporal patterns of IIS and the propagation of the triggering front is associated with the effective hydraulic diffusivity of the subsurface.) Equation 11 suggests that the triggering front propagates with t (Shapiro, 2015;Shapiro & Dinske, 2009b;Shapiro & Müller, 1999;Shapiro et al., 1997). This behavior is occasionally found in IIS data sets and corresponds to a two-dimensional, diffusional process. The upper envelope of the seismic data cluster in Figure 12, though, appears to indicate a different behavior in the propagation of the triggering front. More complex triggering front behavior has been previously explained using nonlinear diffusion in two dimensions (Hummel & Shapiro, 2012Shapiro & Dinske, 2009a, 2009b or by accounting for full poroelastic coupling in a two-dimensional, homogeneous medium (Rozhko, 2010(Rozhko, , 2012Segall & Lu, 2015;Shapiro, 2012). While these physical processes may help to explain the triggering front behavior observed in this data set (and its apparent deviation from the  t r t relationship), we suggest that the triggering front propagation is more appropriately described using a three-dimensional diffusion solution given by Equation 10. The underlying porous medium is in fact three-dimensional and the diffusive propagation of the triggering front according to Equation 10 does not necessarily conform to the  t r t relationship and may therefore provide a better fit with the Soultz-sous-Forêts data set. When fitting Equation 10 to the upper envelope of the seismic data cluster (shown in Figure 12b), we have approximated the fluid source as a point following the approach of Delepine et al. (2004). We have also assumed the fluid viscosity to be 2.0 × 10 −4 Pa⋅s to reflect water at approximately 200 °C and the threshold triggering pressure increment to be 0.04 MPa, which is the minimum value used by Keranen et al. (2014) and well within the range suggested by others (Dempsey & Riffault, 2019;Gischig & Wiemer, 2013;Goebel et al., 2017;Shapiro, 2015).
The triggering front according to Equation 10 fits the upper envelope of the seismic data cluster more accurately than Equation 11. Notice that the rapid spread of seismicity during early times is captured well by Equation 10 while underestimated by Equation 11. The propagation rate of the triggering front at late time appears to be properly estimated by Equation 10 as well. To evaluate the fit of each expression, we calculated both the normalized-root-mean-square-error (NRMSE) and the Nash-Sutcliffe model efficiency coefficient (NSE) for Equations 10 and 11. The reference data set for the location of the triggering front was found by first binning the data into time periods of 0.1 days and then using the 10% of seismic events that occurred farthest away from the injection location within each time period. This set of seismic events are highlighted in Figure 12b in dark gray for reference. The NRMSE for Equations 10 and 11 are 0.05 and 0.09 respectively. The NSE for Equations 10 and 11 are 0.56 and −0.42 respectively. This indicates that Equation 10 is a considerably better fit to the observed data than Equation 11, and may therefore describe the propagation of the triggering front more accurately and produce a more reliable estimate of seismic diffusivity. For the Soultz-sous-Forêts data set, the best estimate of seismic diffusivity was determined to be 4.6 m 2 /s according to the fit of Equation 10 shown in Figure 12b. The potential uncertainty related to this estimate of seismic diffusivity is discussed in Text S5.
Other studies have attempted to characterize the subsurface in the vicinity of the GPK2 well (producing estimates of effective permeability) using various approaches including a laboratory experiment on core samples (Hettkamp et al., 1998), a pre-stimulation slug test (Weidler, 2001), an inverse modeling study using downhole pressure and injection rate data (McClure & Horne, 2011), a circulation test between wells GPK2 and GPK3 (Ledésert & Hébert, 2012), and a tracer experiment between wells GPK2, GPK3, and GPK4 (Vogt et al., 2012). To compare our value of seismic diffusivity to these studies, we first convert the estimates of effective permeability to effective hydraulic diffusivity using Equation 2. For this, we assume that viscosity of the injected water at 200 °C is μ = 2.0 × 10 −4 Pa⋅s, the porosity of the fractured, crystalline rock formation is ϕ = 0.05 (Freeze & Cherry, 1987), and the typical value for the compressibility of water is β f = 4.4 × 10 −10 Pa −1 . For the bulk compressibility of the fractured crystalline rock (β m ), Vogt et al. (2012) employs a value of 10 −8 Pa −1 while Freeze and Cherry (1987) report typical values ranging between 10 −10 Pa −1 and 10 −8 Pa −1 , where values less than 10 −10 Pa −1 pertain to intact crystalline rock. Here, we use a value of 10 −8 Pa −1 when converting to effective hydraulic diffusivity for the estimate associated with Vogt et al. (2012). For the other studies, we estimate a range of effective hydraulic diffusivity using the end member values of bulk compressibility reported by Freeze and Cherry (1987). The resulting estimates of effective hydraulic diffusivity from each study are shown in Table 4. Note that the measurements at relatively smaller scales by Hettkamp et al. (1998) and Weidler (2001) may not be representative of the entire fracture network within the injection formation; therefore, these studies may underestimate the actual effective hydraulic diffusivity at the reservoir scale.
Comparing the value of seismic diffusivity found in this study (i.e., 4.6 m 2 /s) to the values of effective hydraulic diffusivity given in Table 4, it is clear that the seismic diffusivity is one to five orders of magnitude larger than the effective hydraulic diffusivity. This case study thus provides potential empirical evidence suggesting that the seismic diffusivity associated with spatiotemporal patterns of IIS is distinct from the effective hydraulic diffusivity in heterogeneous porous media, such as the highly fractured, crystalline rock formation at the Soultz-sous-Forêts EGS site.

Conclusions and Discussion
In this study, we analyzed the difference between the diffusivities associated with the propagation of seismicity (i.e., the seismic diffusivity) versus mean pressure diffusion (i.e., the effective hydraulic diffusivity) in heterogeneous subsurface formations. Often, the spatiotemporal patterns of IIS are interpreted with r-t plots and the concept of a triggering front, which has been previously shown to propagate in a diffusion-like manner with an associated diffusivity parameter (Segall & Lu, 2015;Shapiro, 2015;Shapiro et al., 1997). Here, we refer to this diffusivity as the "seismic diffusivity" following the terminology of Talwani and Acree (1985). If we assume that the onset of seismicity is induced by a threshold triggering pressure increment (p t ) (Dempsey & Riffault, 2019;Gischig & Wiemer, 2013;Goebel et al., 2017;Keranen et al., 2014;Shapiro, 2015), then the seismic diffusivity is associated with the spatiotemporal evolution of the farthest radial distance to which the threshold triggering pressure increment has propagated (i.e., r max (p = p t )). It is well established that effective hydraulic properties such as the effective hydraulic conductivity and effective hydraulic diffusivity describe the behavior of mean pressure fields (i.e., p avg ) in heterogeneous porous media (Dagan, 1989;Gelhar, 1993;Sanchez-Vila et al., 2006;Zhang, 2001). In a homogeneous porous medium, the seismic diffusivity is undoubtedly equivalent to the effective hydraulic diffusivity. However, due to the rapid preferential propagation of pressure increments through pathways of relatively high hydraulic diffusivity, this equivalence is not expected to hold in heterogeneous subsurface formations, since the propagation of r max (p = p t ) and r(p avg = p t ) are unlikely to coincide.
We have presented numerical simulations of fluid injection and IIS in heterogeneous domains in order to investigate the possible distinction between the seismic and effective hydraulic diffusivity. The fluid flow model (based on uncoupled, linear pressure diffusion) simulated IIS by evaluating the Mohr-Coulomb failure criterion at randomly seeded weak points throughout the domain. We considered two forms of subsurface heterogeneity: spatially correlated, random fields of hydraulic diffusivity (i.e., the SGS domains) and highly fractured, low permeability rock (i.e., the DFNM domain). Results of the SGS simulations show that the location of r max (p = p t ) does indeed exceed r(p avg = p t ) in all cases at all times, with the distance between them increasing with the variability in the hydraulic permeability. Thus, we found that the seismic diffusivity (which describes the propagation of r max (p = p t ) and the spatiotemporal patterns of IIS) is in fact distinct from and greater than the effective hydraulic diffusivity of the SGS domains. After calculating the seismic diffusivity (D s ) and effective hydraulic diffusivity (D eff ) of each SGS domain, we found that the ratio D s /D eff was always above one and systematically increased with the degree of heterogeneity (i.e., with Var(ln κ)). For the largest value of Var(ln κ) considered (i.e., 5.3), the ratio D s /D eff had an ensemble mean of approximately 1.9 and the maximum value encountered in any realization was approximately 3.4. In contrast, results of the injection simulations in the DFNM domain found that the ratio D s /D eff was approximately HAAGENSON AND RAJARAM Note. Ranges shown here account for the broad range in the bulk compressibility of fractured rock (β m ) (Freeze & Cherry, 1987). Overall, the results of this study indicate that the influence of subsurface heterogeneity creates spatiotemporal patterns of IIS that are not accurately described by the effective hydraulic diffusivity of the heterogeneous domain. Rather, the propagation of the triggering front and spatiotemporal patterns of IIS are controlled by the so-call seismic diffusivity (Talwani & Acree, 1985). This result suggests that estimates of hydraulic diffusivity from seismicity-based approaches (Delepine et al., 2004;Hummel & Müller, 2009;Hummel & Shapiro, 2012Rothert & Shapiro, 2003;Segall & Lu, 2015;Shapiro, 2015;Shapiro & Dinske, 2009a, 2009bShapiro & Müller, 1999;Shapiro et al., 1997Shapiro et al., , 2002 likely over-estimate the true effective hydraulic diffusivity of the subsurface. Another critical implication is that modeling the subsurface response to fluid injection operations with homogeneous domains using the effective hydraulic properties of the injection formation (Brown et al., 2017;Catalli et al., 2016;Dempsey & Riffault, 2019;Keranen et al., 2014;Langenbruch & Zoback, 2016;Langenbruch et al., 2018;Pollyea et al., 2019;Riffault et al., 2018) may underestimate the rate of propagation of seismicity and the size of the seismically active region. Alternatively, when employing the seismic diffusivity as a substitute for the effective hydraulic diffusivity, such models would produce the correct spatiotemporal patterns of IIS, yet other hydraulic processes may be inaccurately represented (such as inter-well connectivity or reservoir pressurization) (Birdsell et al., 2018;Haagenson et al., 2018).
Even if seismic diffusivity is not an accurate estimate of the effective hydraulic diffusivity in heterogeneous porous media, it may still be helpful by providing a simple descriptor of the rate at which seismicity appears to spread at a particular location. Since both the seismic diffusivity and the effective hydraulic diffusivity of subsurface formations are influenced by heterogeneity, estimation of the statistical properties of the underlying heterogeneity based on estimation of the seismic diffusivity may provide an indirect approach to estimating the effective hydraulic diffusivity. Conversely, hydraulic characterization at a particular site may indirectly provide estimates of the seismic diffusivity, allowing operators to evaluate the potential for the rapid propagation of seismicity.
Our approach involves approximations, which implies certain limitations. Foremost, the subsurface stress state is considered static during the simulation and only impacts the calculation of the Mohr-Coulomb failure criterion. This approach is similar to numerous previous studies that also neglect potential mechanical effects on the behavior of IIS (Brown et al., 2017;Dempsey & Riffault, 2019;Hummel & Shapiro, 2013;Keranen et al., 2014;Langenbruch et al., 2018;Nakai et al., 2017;Rothert & Shapiro, 2003;Shapiro, 2015;Shapiro & Dinske, 2009b;Shapiro et al., 1997Shapiro et al., , 2002Talwani & Acree, 1985). Recently, there is a growing body of research investigating the influence of mechanical coupling on the behavior of IIS through both poroelastic stressing (Chang & Segall, 2016;Jha & Juanes, 2014;Rutqvist et al., 2013;Segall & Lu, 2015;Zhai & Shirzaei, 2018;Zhai et al., 2019) and static stress transfer following a seismic event (Catalli et al., 2016;Schoenball et al., 2012), which will certainly lead to improved insights on the physical nature of IIS phenomenon. In addition, we have neglected the dilation of fractures due to the change in normal effective stress (Bandis et al., 1983) or due to shear failure along a fracture (Rong et al., 2016;Ye & Ghassemi, 2018). We fully expect that including fracture dilation would only hasten the propagation of the triggering front and thus enhance the distinction between seismic diffusivity and the effective hydraulic diffusivity in formations of fractured rock. The influence of mechanical coupling is more difficult to predict, as changes to the stress state could be either stabilizing or destabilizing depending on the orientation of any given fracture. Future studies on this topic could investigate the potential influence of mechanical coupling or fracture dilation on the spatiotemporal behavior of IIS in three-dimensional domains of fractured rock. It may also be insightful to further explore the critical implications of the study's findings-particularly the potential for connecting the spatiotemporal patterns of IIS and estimates of the seismic diffusivity to various properties of the underlying fracture network. This could provide a highly beneficial tool for subsurface characterization or predicting the seismic response at potential fluid injection sites. Fluid mass source or sink term per unit area of fracture Q i Total fluid flux in i th direction of numerical permeameter test r Radial distance from injection location r t Radial distance from injection location to triggering front at time t r max (p = p t ) Farthest radial distance to which threshold triggering pressure increment propagates r(p avg = p t ) Radial distance at which directionally averaged pressure increment equals threshold triggering pressure increment t e Duration of fluid injection t t Triggering front arrival time at given radial distance from injection location (

Data Availability Statement
The seismic data set used in the case study of Soultz-sous-Forêts is freely available through the EPOS-IP Anthropogenic Hazards data repository (Leptokaropoulos et al., 2019).