High Injection Rates Counteract Formation of Far‐Reaching Fluid Migration Pathways at The Geysers Geothermal Field

Deep underground water injections induce seismicity. When the seismic fractures coalesce into far‐reaching pathways for fluid migration, the migrating fluid may reach preexisting faults and, by decreasing fault strength, can trigger major seismic events. We assume that the potential for building such pathways depends on closeness of hypocenters, similarity of fracture planes orientations, and closeness of radii taking off from the injection point, on which events locate. We define this potential as the average distance between seismic events in the space of parameters quantifying the above conditions. We show that in the studied case from The Geysers geothermal field, this potential is highly correlated with injection rate. When the overall level of injection rate is high, the higher the injection rate is, the more the potential for building far‐reaching pathways for fluid migration is reduced.


Introduction
Deep underground water injections induce seismicity. This seismic fracturing of rocks is desirable as it increases the surface on which heat exchange takes place in Enhanced Geothermal Systems. However, the seismic fractures may also coalesce into undesired pathways for fluid migration. These are such pathways that enable the fluids to reach preexisting faults. In result, by decreasing fault strength, the fluids can trigger ruptures and produce major seismic events. The further the migration pathways extend from the injection point, the more probable are these unwanted effects. It is therefore of paramount importance to recognize under which injection conditions an induced seismic process can produce such pathways (e.g., Davies et al., 2013;Ellsworth, 2013;Majer et al., 2012;Zhang & Yang, 2015).
Numerous studies indicate that the geometry of fractures and the structure of fracture networks are the main factors controlling fluid flow and the fluid transport characteristics of rocks (e.g., Hope et al., 2015;Lee et al., 1990;Long & Billaus, 1987;Schwartz et al., 1983;Snow, 1965). The development of the fracturing process has been investigated in a number of laboratory experiments (e.g., Ko & Kemeny, 2011;Lockner et al., 1992;Stanchits et al., 2006). Based on these experiments, many models have been developed to simulate the geometry of fractures and the topology of fracture networks (e.g., Hope et al., 2015;Lee et al., 1990;Long & Billaus, 1987). However, although seismic field data provide direct insight into the development of fracture networks at the crustal scale, there are few studies focused on this topic (e.g., Chorozoglou et al., 2018;Kagan, 1992Kagan, , 2000 Here we study seismic and injection data from a part of The Geysers geothermal field to determine a relationship between the injection conditions and the potential of injection-induced seismicity to build far-reaching pathways for fluid migration. We formulate three conditions which we expect to play a role in linking fractures and building such pathways: closeness of hypocenters; similarity of fracture planes orientations; and closeness of radii, which begin at the open hole section of the injection well and on which events occur. We assume that in the same injection conditions and for the same orientation of the line connecting hypocenters of two events with respect to the orientation of regional stress field, the probability for these events to link is higher when they are closer to each other than when they are farther from each other. We assume that, for the same stress and injection conditions and the same distance between hypocenters, when the fault planes of two events are parallel and they are parallel to the line connecting hypocenters, the probability for these events to link is higher than this probability for other mutual orientations of fault planes. Moreover, when they have linked, they more likely extend farther than the linked fractures with other fault plane orientations. We assume that linked fractures located along the straight line beginning at the injection point reach farther from this point than such fractures located in another way. Consequently, seismic events are represented by eight parameters: three hypocentral coordinates, three independent angles determining orientations of the T and P axes of the double-couple focal mechanisms, and two angular coordinates of hypocenters in the spherical system beginning at the open hole of injection well. To achieve the same scaling of these parameters, we transform them to equivalent dimensions (EDs) (Lasocki, 2014; Supporting Information Text S1). The average distance between the events in the eightdimensional space of the aforementioned parameters, called the degree of disordering of sources, ZZ, expresses to which extent the above three conditions have been fulfilled. The chance for the seismic events with small value of ZZ that they link and reach far is higher than in other cases. In this way, ZZ quantifies the potential for building far-reaching fluid migration pathways.
We show that ZZ is highly significantly correlated with injection rate. ZZ took the highest values when the injection rates were the highest, which indicates that high injection rates counteracted the formation of farreaching fluid migration pathways.

Data and Methods
The Geysers geothermal field is located in California, USA. Geothermal operations, which began there in 1960s, and which presently use EGS technology, have induced hundreds of thousands of seismic events. We studied injection and seismic data from an isolated area of 2 km × 2 km in the NW part of The Geysers, between 10 December 2007 and 23 August 2014. The basis for the seismic data set was an improved catalog that contained 1,252 events (Kwiatek et al., 2015;Martínez-Garzón et al., 2014, 2016. This catalog provided all the event parameters, necessary for the calculation of the degree of disordering of sources, ZZ: hypocenter locations with an accuracy of 50 m and focal mechanisms with an accuracy of 20°/5°/10°for strike/dip/rake. The focal mechanisms were used to recover the trend and plunge angles of T and P axes. Moment magnitude was used, and the completeness level was Mc = 1.4 (Kwiatek et al., 2015). The strongest event was Mw3.2. Hypocenter locations and magnitude distributions are provided in Figures S1-S5.
In the study period, the injections in the study area were carried out into two wells: Prati9 and Prati29. There were three phases of the injection activity: 1. Phase F1 from 10 December 2007 to 10 April 2010, in which only Prati9 was operational; 2. Phase F2 from 11 April 2010 to 21 June 2013 with simultaneous injections into both wells; and 3. Phase F3 from 11 June 2013 to the end of the study period, in which only Prati9 was operational.
The injection data consisted of the daily injection volumes into Prati9 and Prati29.
Most of the 1,252 seismic events clustered around the open hole of Prati9 well (see Figures S1-S4). We intended to analyze the relationship between injection and the degree of disordering of sources, ZZ, within the whole time period of available data, in which period only Prati9 well was constantly operational. Therefore, we studied these events located closer to Prati9 well, expecting that their spatial relation to Prati9 resulted from their physical relation to the injection activity in Prati9. As the selection criterion, we used the hypocentral distance from the open hole of Prati9 to be no more than 600 m; 1,121 events that fulfilled this criterion are called cluster A.
The degree of disordering of sources, ZZ, used to quantify the potential for building the far-reaching fluid migration pathways was the average distance between the seismic events in the eight-dimensional parameter space {x 1 , x 2 , x 3 , plu_X 1 , plu_X 2 , tre_X 1 , Θ, φ}. x 1 , 2 , 3 were hypocenter coordinates. plu_X 1 , tre_X 1 were the plunge and trend angles of T axis. plu_X 2 was the plunge angle of P axis. Θ, φ were the polar and azimuthal angles of hypocenter in the spherical system of coordinates beginning at the open hole of Prati9 well. Θ = 0 for the vertical direction, and φ = 0 for the N direction. These parameters were not comparable; therefore, we first transformed them to ED. The transformation to ED is a technique based on a probabilistic equivalence of the parameters that scale differently (Lasocki, 2014; Text S1). All transformed parameters are uniformly distributed in [0, 1], and the distance between any two objects is the Euclidean metric. Further on the symbols x 1 , x 2 , x 3 , tre_X 1 , tre_X 2 , plu_X 1 , plu_X 2 , Θ, φ denote the EDs of the original source parameters.
The closeness of hypocenters was parameterized by the absolute differences between hypocenter The trend angle and the polar angle take values in [0°, 180°], which is [0, 1] in ED, and the azimuthal angle takes values in [0°, 360°], which is also [0, 1] in ED. The shortest distances between these angles in ED were therefore The multipliers in front of the opening braces in equations (1)-(3) were inserted so that the differences in all parameters scaled in [0, 1].
For a collection of n seismic sources, the degree of disordering of sources, ZZ, reads: where where Δ r (i,j) is the distance between hypocenters of events i and j; Δ M (i,j) is the distance between focal mechanisms of these two events; and Δ ϕ (i,j) is the distance between the directions of radii from the Prati9 open hole, those on which the hypocenters of these two events locate.
We carried out our analyses separately in the three injection phases. Out of 1,121 studied events, 248 events occurred in the injection phase F1, 702 events in phase F2, and 171 events in phase F3. For every injection phase, we calculated ZZ for 50-event window sliding by 10 events because the number of events per phase was not high enough to use non-overlapping windows. We obtained 21 ZZ values for F1, 66 ZZ values for F2, and 13 ZZ values for F3.
Next, we calculated the average injection rates, IN, during time windows covering the periods of the 50-event sliding windows, respectively. The calculations were performed for the injection rate periods exactly matching the periods of respective 50-event windows and for the injection rate periods from 1 to 21 days preceding the periods of event windows. The latter part of this analysis was meant to study delayed reactions of seismicity.
Finally, we studied the correlation between ZZ and IN. It comes from equation (4a) that ZZ is composed of three components representing our three conditions determining the potential of injection-induced seismicity for building far-reaching pathways for fluid migration. In order to recognize contributions of these components to the correlations between ZZ and IN, we studied also the correlations between Δ r ,Δ M ,Δ ϕ and IN, respectively.
All correlation analyses were preceded by the Jarque-Bera test applied to check normality of the distributions of used variables. If the normality hypothesis was not rejected, we used Pearson correlation; otherwise, we used Spearman correlation.
We also tested differences between the values of location parameters by phase of IN, ZZ, seismic activity rate, and magnitude, respectively. When the normality hypothesis for a parameter was not rejected, we used the Student's t test for means; otherwise, we applied the Mann-Whitney U test for medians. All statistical inferences were performed under the significance level α = 0.05.

Results and Discussion
Some descriptive statistics by injection phase of IN, ZZ, seismic activity rate, and magnitude of events are presented in Table 1 together with comparisons of their location parameters (means or medians) among phases. The medians of injection rates into Prati9 well were decreasing statistically significantly with injection phase. . It is seen that in the first two injection phases, F1 and F2, ZZ correlated positively with IN. Also, the amplitudes of the ZZ changes agreed well with the amplitudes of the average injection rate changes. In phase F2, this agreement related to the summed injection into both wells (blue curve) rather than into Prati9 well alone (black curve), even though the analyzed seismic events were geometrically linked to Prati9 (Figures S1-S4). correlations with Δ ϕ in phase F2 because Δ ϕ series in this phase was not normally distributed. In these two cases Spearman's rank correlation was used.
The results in Table 2 confirm the relationships evidenced in Figure 1. In F1 and F2 the correlation between IN and ZZ was significant, positive. In F2 this correlation was highly significant, irrespective of whether the IN referred to injections into Prati9 or to the summed injections into Prati9 and Prati29 wells. The ZZ versus IN(9) scatterplot is presented in Figure S6.
The results of correlation analysis for delayed ZZ with respect to IN are presented in Tables S1-S4. In F1 the correlation coefficient slowly decreases with increasing lag. The correlation coefficient in F2 slightly increased when ZZ was delayed with respect to IN; however, the difference between the correlation coefficient for zero lag and for 13 days lag was only 0.022 and was surely insignificant.
The previous studies of the same data have indicated a positive correlation between seismicity rates and injection rates (Leptokaropoulos et al., 2018). Hence, different time periods corresponded here to the 50event windows; the time periods at higher injection rates were shorter. We checked whether the positive ZZ versus IN correlations were not due to these differences between time periods. We divided phase F2 into non-overlapping windows of constant, 68-day duration. The correlation between IN and ZZ both calculated for these windows of constant time period was also positive (0.87) and significant (p = 0.0097).
In F1 only Δ r out of three components of ZZ significantly and positively correlated with IN. Hence, in this injection phase, the positive ZZ-IN correlation resulted from that that higher injection rates were increasing distances between the sources.  In F2 all three distances, Δ r ,Δ M ,Δ ϕ , were highly positively correlated with IN; thus, they all significantly contributed to the correlation ZZ-IN. Higher injection rates led to an increase of the distances between hypocenters, to a greater variety of P and T axes directions, and to a greater angular dispersion of the hypocenters in relation to the open hole of Prati9 well.
The significant correlation Δ r -IN in F1 and F2 cannot be attributed just to moving away sources from the well opening of Prati9 when the injection rates were higher. The correlation between the average distance from the Prati9 opening and IN was significant in F1 (corr. coefficient 0.47, p = 0.03), though much weaker that the Δ r -IN correlation, and in F2 it was insignificant (corr. coefficient 0.14, p = 0.28).
Our target was studying the relationship between the degree of disordering of seismic sources, ZZ, and the injection rate, IN. We used the transformation to ED because the components of ZZ, the distances between P (and T) axes of events and the distances between radii on which events located, are not Euclidean. However, the metrics of the distance between hypocenter locations is Euclidean. Thus, it was possible to compare the shown here results of correlation between Δ r and IN in F1 and F2, based on the transformation to ED, with such results based on more often used inter-event distance analyses. As alternatives to Δ r , we used correlation dimension of distances between hypocenters, D 2 , and summarized squared distances in the original (x,y, z) space, d2. D 2 turned out to be uncorrelated with IN, probably because of difficulties with its estimation. The correlation between d2 and IN was even stronger than between Δ r and IN. We think that this was due to the fact that any transformation, here the transformation to ED, causes some loss of information. Nevertheless, the same significant correlation between d2 and IN as the correlation between Δ r and IN positively validated the performance of transformation to ED. The results of this part of analysis are provided and discussed in Text S2.
The significant impact of fluids on the stress field and the faulting regime is known from many studies (e.g., Bachmann et al., 2012;Hardebeck & Hauksson, 1999;Segall & Fitzgerald, 1998). The analyses of seismic events from the NW part of The Geysers geothermal field, presented in Martínez-Garzón et al. (2013, 2016 and Kwiatek et al. (2015), clearly show a large variability of focal mechanisms. Explaining the observed stress tensor perturbation, the cited authors consider the fact that in addition to fracture reactivation, massive fluid injection in EGS systems results in hydro-fracturing. It is then possible that during the time periods of higher injection rates, new small fractures were created, and they could also perturb the stress field in the observed way (Martínez-Garzón et al., 2016). A number of events either occurred on severely misoriented faults (low instability coefficient, e.g., Vavryčuk, 2011) or slipped in a different orientation than that predicted from the stress field. These events mostly occurred during periods of high injection rates indicating that faults not optimally oriented to the stress field require larger pore pressures to become activated.
A way in which increased injection rates may increase the degree of disordering of sources (ZZ) is shown schematically in Figure 2. An increase of pore pressure resulting from the increased injection rate broadens the range of possible orientations and locations for new shear fractures.
Studying fracture network development on the same data as here, Orlecka-Sikora et al. (2019) have shown that the connectivity of fractures induced by fluid injection is lower for higher injection rates and vice versa. Thus, the connectivity was responding to the injection rate changes in the same way as ZZ. This allows for inferring that ZZ indeed represents the potential for fluid migration. In the cited work the connectivity was estimated by the connectivity coefficient, C, defined by the ratio of the observed number of intersections of fractures in a fracture network to the number of all possible intersections in this network. These all possible intersections were evaluated as 0.5f(f − 1), where f was the number of fractures building the fracture network (Albert & Barabasi, 2002;Watts & Strogatz, 1998).
Prati29 well was operational only in phase F2. In this phase the amplitudes of ZZ changes agreed with the amplitudes of changes of the summarized injections into Prati9 and Prati29. Furthermore, the mean level of ZZ and the mean seismic activity were the highest in this phase and significantly higher than in phases F1 and F3. In F2 only the level of total injection into Prati9 and Prati29 was the highest, and the mean level of injection rate into the Prati9 well was significantly smaller than in phase F1. All these indicate the important role of the injections into Prati29 for the generation of seismic events within the studied cluster A despite the Prati29 well location being outside cluster A.
The seismic events generated at higher injection rates, that is, the ones more disordered (ZZ was higher), did not have greater magnitudes. On the contrary, in phase F1 the mean magnitude in windows correlated negatively with ZZ (corr. coefficient −0.61, p = 0.003). In F2 no magnitude-ZZ correlation was ascertained.
In phase F3, in which the overall level of injection rate was the lowest among injection phases, the correlation IN-ZZ was significant, negative. This correlation was achieved only jointly by the three components of ZZ, Δ r ,Δ M ,Δ ϕ , because neither of them significantly correlated with IN. The negative correlation coefficient in F3 increased when ZZ was delayed with respect to IN, and for 4 and more days lag, it became statistically not significant (Table S4). This reversal of correlation results may be connected with another fracture mechanism below a certain level of injection rate. We discuss this possibility in Text S3. However, the data series in F3 was composed of only 13 points; therefore, the correlation (p ≈ 0.03) might be spurious.

Conclusions
Studying the actual field data, we found exceptionally high and immediate correlation between the injection rate and the seismicity parameter-the degree of disordering of sources, ZZ. Also the amplitudes of ZZ changes agreed well with the amplitudes of average injection rate changes. ZZ is defined solely on parameters of seismic sources. It describes how much seismic fractures are dispersed in terms of distances between their hypocenters, mutual orientations of their fracture planes, and angular dispersion of their hypocenters. In the next work we will modify ZZ to account also for source sizes.
We interpret ZZ as a measure of the potential of seismic fractures for building far-reaching pathways for fluid migration. The logic of the three conditions, on which ZZ is based, supports this interpretation. However, ZZ does not represent all possibilities for building such a fracture system neither it is adequate for other possibilities of fracture network development (e.g., not for a percolating fluid pathway).
In the studied case from The Geysers geothermal field, the optimal conditions to avoid such ordering of seismic fractures that enable linking them into longer pathways, extending farther from the injection point, were met for high injection rates. The higher the injection rate was, the more disordered the seismic fractures were generated; that is, the chances to build longer pathways for undesired fluid migration decreased. High injection rates caused an increase in seismic activity; nevertheless, the median level of seismic event magnitudes remained unaffected.
The above conclusion, if confirmed in other cases of injection-induced seismicity, would help to understand this type of seismicity and related hazards. Martínez-Garzón et al. (2018) concluded that the used here data set from the NW region of The Geysers well represents the broader seismic processes at The Geysers field. Thus, our results may be also valid for all seismic processes at The Geysers. However, for further generalizations studies of other injection-induced seismicity cases are required.