Evidence for subcritical rupture of injection-induced earthquakes

Seismicity induced by geo-engineering operations may be hazardous for people, infrastructure and the environment. The crucial information for assessing induced seismic hazards and related risks is knowledge of the time-dependent strength of rocks and the deformation due to fluid injection. Our studies of seismic and injection data from a geothermal field indicate that pressurized injections lead to rock fracturing at stress levels below the rock toughness, i.e., subcritical fracture growth. We provide a relation between the rate of this subcritical fracture growth and the injection rate. Based on this relation, we estimate the maximum subcritical magnitude. We hypothesize that subcritical fracture growth may be controlled by the amount of stress asymmetry, i.e., the relative values of the principal stresses. We discuss the conditions under which the subcritical fracturing regime can transform to a critical state and critical rupture may occur. We present the possibility of using these results in the operational reservoir to manage seismic hazards.

www.nature.com/scientificreports www.nature.com/scientificreports/ fracturing process, we propose to modify the formula of Galis et al. 8 by incorporating a parameter describing the rate of the fracturing process into the formula. We also discuss the stress state under which the fracture networks grow subcritically. We find that the subcritical rupture regime can be controlled by the stress parameter R ( )/( ), , are the principal stresses. We find that if 2 σ is close to σ 1 or 3 σ , the rate of SFG is statistically significantly lower than for the conditions when the principal stresses are more symmetrical, i.e., σ 2 is centred between σ 1 and 3 σ . Finally, we analyse conditions under which the subcritical rupture regime can experience a transition to a critical one and critical rupture may occur.

Data
We investigate the development of fractures using seismic data recorded in the northwestern part of TG geothermal field in California, USA, related to large-scale, long-term fluid injection into two wells, Prati-9 and Prati-29. TG is located in a series of complex, intensely deformed and faulted metamorphic Franciscan greywackes (e.g., 9 ). The seismicity in this region results from the thermoelastic and poroelastic effects that influence the local stress field in the vicinity of the injection wells (e.g., [9][10][11][12] ). The regional stress regime in the TG field is extensional. Based on the fault plane solutions of seismicity recorded in this area, a NNE orientation of the maximum horizontal stress is found, indicating a combined normal and strike-slip faulting regime, and the strike-slip component increases with reservoir depth 13 . The catalogue of 1254 seismic events and injection rate data is available on the induced seismicity European Plate Observing System platform of the EPOS Thematic Core Service Anthropogenic Hazards (IS-EPOS, https://tcs.ah-epos.eu). Seismic source parameters of these events were calculated by Martínez-Garzón et al. 9 and Kwiatek et al. 10 For the purpose of the following analysis, the stress drops are calculated based on seismograms using seismic moments and corner frequencies following 14 (Materials and methods).
We consider two of the injection cycles described in 9 . We select these cycles after 9 since they frame the most intensive fluid injection into the reservoir, resulting in a large number of seismic events. In total, 509 seismic events occurred in the period between April 2008 and October 2011 10 . The cycles include peaks of fluid injection during two different periods, i.e., when only the Prati-9 well is active and when both wells Prati-9 and Prati-29 are operating. Moreover, detailed technological activity data during these periods are available in the literature (e.g., 9 ). Each cycle is divided into three stages: preceding, during, and following peak injection. The selected cycles exhibiting these stages are presented in Fig. 1, along with the daily injection data.
Our considerations of SFG refer to system-sized ruptures; thus, we use the fracture networks identified in TG by Orlecka-Sikora et al. 15 . Because the reservoir rocks in TG are highly fractured, the authors assumed that its fractures propagate along pre-existing fractures, which has also been suggested to be the dominant mode of failure at the crustal level 2 . Furthermore, the authors used seismic source parameters to identify similarly oriented fractures according to criteria based on the fault plane orientations of fractures and the locations of hypocentres with respect to the injection well and the regional principal stress orientation. These fractures, referred to here as a fracture network (FN), are assumed to represent the damage in the vicinity of the same system-sized ruptures. To identify fracture networks, these authors used hierarchical clustering. This approach resulted in the identification of 13 FN in the time period we are focused on here 15 .

Results
In subcritical hydraulically driven fracture growth, a fracture grows due to either an increase in fluid pressure or a decrease in normal stress (e.g., 16,17 ). The long-term loading of pore fluid stresses and thermal stresses can weaken a rock's resistance to fracture. These mechanisms have been shown to be responsible for the seismicity observed in TG (e.g., 9,10,18 ). Previous studies of seismic moment tensors in TG revealed a mixed-mode fracture mechanism [19][20][21][22] . Laboratory experiments have determined that the values of the parameters for mode II and mode III SFG are similar to the corresponding values of mode I SFG regardless of the loading configuration or specimen geometry 4,23 . Based on these results, we assume that the constitutive equations describing subcritical tensile crack growth hold for all three fundamental modes of loading (e.g., 2,4,23 ). The most widely used model is Charles' power law 24 , which describes the crack tip velocity in a subcritical regime 2 . Charles' law can be considered in the following form, which describes the rate at which a fracture grows: a where l is the fracture length or diameter; a is the growth exponent, which is related to the stress corrosion index n by = a n/2, and C is a parameter that depends on the stress state. The acceleration of the fracture length is predicted by solving Eq. (1): where l 0 is the initial fracture length at time = t 0 and t f is the failure time 25 . In the case of a constant stress, instability can develop if > n 2. Therefore, if the average increase in fracture unit length over time is high, it would lead to sudden unstable fracture propagation, resulting in a runaway earthquake.
The fracture length is proportional to the cube root of the seismic moment, M 0 26 . The relation between M , 0 and the seismic source dimension, r , 0 for the analysed dataset was provided by Kwiatek et al. 10 The moment magnitude is calculated from M 0 using the equation of Hanks and Kanamori 27 : . We focus on the fracture network growth rate. The criteria applied to identify fracture networks by Orlecka-Sikora et al. 15 www.nature.com/scientificreports www.nature.com/scientificreports/ allow us to track the coseismic increases in fractures that occur during seismic events. In our considerations, rupturing due to seismic events builds a fracture network. Such an approach is in line with the observations of injection-induced seismicity and the results of laboratory experiments (e.g., 28,29 ). However, in the case of geothermal reservoirs, there is evidence for the importance of aseismic growth processes (e.g., [30][31][32][33]. We assume that the ratio of seismic to aseismic growth is constant under certain stress conditions; thus, we take the total seismic length as a proxy for the total seismic and aseismic length. Hence, the complexity of fracture growth directly reflects the complexity of the developing network structure. We consider the growth rate for every FN identified by 15 during each injection stage in both cycles. To reduce the data scatter, especially that in the injection rate loading data, when considering dl dt , we use the observed average increase in FN length per unit time during a particular stage of injection as follows: where i n 1, = … is the number of observed increases in the FN in the ∆t time period due to constant loading, which is assumed using the mean injection rate (mIR) in ∆t. The FN propagates with a mean velocity of v mean in t ∆ , and l is the total length of the FN at time t. The SFG parameters a and C are determined based on the slope and intercept of the linear regression through all v log( ) mean t ,∆ values versus l log ( ) values. evidence for SfG. We observe two types of patterns for − v l log( ) log( ) dependency; the most common one is the linear pattern, and the second pattern is close to the experimentally derived characteristic K v − diagram for tensile SFG, which has three distinct regions of − K v dependency 2 , where K is the stress intensity factor (Fig. 2c). Generally, the stress intensity factor at the crack tip is proportional to the applied stress and the square root of the fracture length (Eq. (3)) 4,34,35 . The complete K v − diagram for tensile SFG derived from studies of glass is assumed to hold for all three fundamental modes of crack displacement, although there is little evidence to support this assumption 2 . The behaviour in region 1 (the area marked as 1 in Fig. 2c) is assumed to be controlled by the rate of stress corrosion reactions at the crack tips. Region 2 (marked as 2 in Fig. 2c) is controlled by the rate of transport of reactive species to the crack tips. In region 3 (marked as 3 in Fig. 2c), crack growth is www.nature.com/scientificreports www.nature.com/scientificreports/ mainly controlled by mechanical rupture 2 . Most experimental data obtained from studies of tensile SFG in rocks appear in region 1 or region 3 on the − K v diagram. Region 2, which corresponds to the diffusion process, is very rarely observed in rocks 2 .
Of the 66 fracture network growth periods associated with the stages of a particular injection cycle, we can estimate the values of the growth exponent a and the parameter C for 33 cases. Information about the uncertainties of these estimates is provided in the Materials and methods section. Although the standard errors of estimation (SEEs) of the a and C parameters are very low for most of the injection stages, we observe variations in the a and C parameters related to the injection rate (Fig. 3). The results reveal a cyclic-like pattern in the relation between the growth exponent, a, and the mean injection rate, mIR, in the stage. In the first range of mIR values (Pattern 1 in Fig. 3), which reach up to approximately ⋅ 7 10 m /day 3 3 , the relation between the injection rate and growth exponent is = . + . ⋅ a m IR 0 15 0 1 . The Spearman correlation coefficient between the injection rate and growth exponent is 0.9 and is statistically significant. Then, at injection rates ranging from -. . ⋅ 4 5 7 5 10 m /day 3 3 (Pattern 2 in Fig. 3), the relation is repeated with higher scatter, with a statistically significant correlation coefficient of 0.7. When mIR exceeds 8 0 10 m /day 3 3 . ⋅ (Pattern 3 in Fig. 3), we observe a relation in mIR versus a with a statistically insignificant correlation coefficient of 0.3. When mIR is standardized within Patterns 1 to 3, the relation of mIR and a is statistically significant with a correlation coefficient of 0.5; the value of a is equal to approximately .  (c) Schematic relation between the stress intensity factor and crack velocity for subcritical tensile crack growth; K IC is the critical stress intensity factor, K 0 is the stress corrosion crack growth limit, and 1-3 represent different behaviours of SFG. www.nature.com/scientificreports www.nature.com/scientificreports/ Seismic potential of SfG in tG -the maximum subcritical magnitude. Since injection-induced seismicity in TG provides evidence for SFG, we analyse its seismic potential from the model of the maximum magnitude of self-arrested rupture, M w max arr − , proposed by Galis et al. 8 We apply this model to each step of FN growth due to an increase in the injection rate. According to 8 , the maximum arrested seismic moment, M max arr 0 − , and the corresponding − M w max arr can be estimated by approximating the additional stress drop induced by the perturbation of pore pressure as a point load superimposed on the background stress drop (i.e., the stress drop before injection), ∆σ 0 . This approximation yields the dependence of earthquake size on the injection volume, V ∆ , with an exponent of 3/2 8 : where κ is the bulk modulus of the reservoir rock, µ d is the dynamic friction coefficient of the fault and h is the reservoir thickness. In this model, we use parameter values that are plausible for TG: the bulk modulus of the fracture-dominated reservoir rocks k Pa 3 6 10 9 = . ⋅ , Poisson's ratio 0 25 ν = . , = h m 1500 , and 0 1 d µ = . 9,36 . We assume after Galis et al. 8 that earthquake stress drops reflect the background stress drop, i.e., the stress drop before the following injection loading, ∆σ 0 . Hence, we apply the mean ∆σ 0 in each injection stage, following the stage for which − M w max arr is calculated. The volume of injected water is calculated based on mIR and the duration of the stage in a particular FN. Generally, the temporal evolution of the maximum observed magnitude follows the behaviour predicted by the model of Galis et al. 8 (Fig. 4). However, the difference between the upper limit predicted by 8 and the observed maximum magnitude increases with time. During the first cycle of water injection, the observed maximum magnitudes for the particular stages of fracture networks correspond well to the estimated M w max arr − . However, during cycle 2, the observed maximum magnitudes are much lower than the arrested model estimates. This observation is not unique; the same trend can be observed for the recorded maximum magnitudes during a 6.1-km-deep geothermal stimulation in Finland 37 . Previous studies of the deformations and seismicity in TG showed that the bulk modulus varies in a wide range (e.g., 36 ). The variability of the bulk modulus for the reservoir at TG can be associated with the fracture density and is lower for more fractured rocks (e.g., 36 ) and higher when pressure increases (e.g., 38 ). Since the SFG exponent provides insights into the fracturing process, we propose to modify the formula of Galis et al. 8 by incorporating the a parameter to estimate the effective bulk modulus for TG: the higher the a value is, the lower the k modulus. Based on the well-confirmed functional dependence of the rock bulk modulus on crack density (e.g., 38 ), we use a simple coefficient = − ⋅ . ⋅ k a P a (1 ) 3 6 10 9 to obtain the effective value of the bulk modulus for the rocks experiencing a particular FN. Because the relation − mIR a is shifted with mIR for Patterns 2 and 3 ( Fig. 3), which is linked to the pressure changes because the fracture density increases, the − a (1 ) coefficient is reduced by (1-shift). Figure 4 presents the maximum magnitudes observed during a particular stage of water injection together with the estimated values of − M w max arr from the seismic moment according to Eq. (4) and with the modified bulk modulus. The differences between the observed maximum magnitudes and upper limit derived from Galis et al., 8  the instability. A fracture growing subcritically can reach instability and become the critical fracture propagating dynamically. The instability requires a stress corrosion index > n 2. In our studies, the a parameter is equal to the stress corrosion index n divided by 2. Moreover, our modified approach of tracking changes in the www.nature.com/scientificreports www.nature.com/scientificreports/ mean velocity of fracture network growth instead of using absolute velocity data (Eq. (3)) influences the a value (Materials and methods) by approximately 30%. After recalculating a values to the stress corrosion index n and correcting for the mean velocity approach, we find some stages of fracture networks growing under > n 2 (Fig. 5). Figure 5 presents the instabilities together with the strongest earthquakes recorded in the considered dataset. The occurrences of the strongest earthquakes in the stages are well correlated with the occurrences of instabilities. Fracture network 1 during the 3 rd stage of the 2 nd injection cycle experienced growth acceleration, but the strong event with M 2 5 w . occurred 3 months later, after the analysed time frame of the 2 nd cycle. The U Mann-Whitney test, which tests whether two independent samples were selected from populations having the same distribution, confirms that the fracture networks with approximated values of parameter n > 1.6 produce higher maximum observed magnitudes than the FN stages with lower n, with the Z statistics value of the test equal to 2.6 for a sample size of 15 for both groups, at a significance level p 0 009 = . .

Stress changes influencing SFG in TG.
Thermo-and poro-elastic stresses are mainly responsible for the seismicity observed in TG. However, both mechanisms can act at the same time to decelerate or accelerate fracture growth 1 . These mechanisms can thus increase the crack density, leading to an increase in the growth parameter a due to crack-linking processes that facilitate fracture extension. Alternatively, they can also lead to a decrease in microcrack density, thus inhibiting SFG. To infer the ability of rocks to develop microcracks, we focus on stress symmetry breaking. The degree of stress symmetry breaking yields information about how far the value of the intermediate stress ( 2 σ ) is from the midpoint between the values of the maximum principal stress ( 1 σ ) and the minimum principal stress (σ 3 ). The results of laboratory experimental studies show that the effect of dilatancy is influenced by the intermediate principal stress (e.g., 39,40 ). If the intermediate stress is close to 1 σ or σ 3 , the dilatancy increases. However, if the principal stresses approach symmetry, i.e., 2 σ is centred between s1 and 3 σ , the rock is strengthened, leading to a significant weakening in the effect of dilatancy. For this purpose, we calculate the relative stress magnitude = following the stress inversion methodology used in the STRESSINVERSE software package 41 . STRESSINVERSE calculates the stress orientation based on focal mechanisms (strike/dip direction/dip angle). The inversion is based on Michael's method 42,43 , in which an instability criterion proposed by Lund and Slunga 44 is incorporated. We use 1000 noise samplings of the input focal mechanisms to estimate the 95% confidence intervals of the relative stress magnitude. We perform a temporal analysis of the stress field changes and assess the impact of depth on the stress inversion results. For the temporal analysis, we apply the stress inversion to the entire dataset, i.e., without dividing the data based on the fracture networks due to the limited number of seismic events in the stages of the FNs. Then, we perform a stress inversion of the focal mechanisms from the moving windows of 50 earthquakes using a step size of 1 event to detect small variations in the relative stress magnitude. The number of events in the moving window is selected to achieve balance via a trade-off between the discrimination of different injection rates and the requirement of a certain variety of focal mechanisms. The rock strength is expected to be highest at = . R 0 5 and lowest at = R 0 and R 1 = . Hence, we analyse the absolute value of this difference = − . AR R 0 5 as a measure of the amount of stress symmetry breaking, with AR 0 = representing the maximum rock strength. Although the R parameter is strongly uncertain, its relative changes are plotted versus mIR (Fig. 6). We observe both rock strengthening and rock weakening when decreasing and increasing the injection rate.
Since the AR value is calculated for the entire dataset in a moving window in time, directly assessing the relation between the a parameter and AR for the FN stages is difficult. To enable this determination, we select in each window the most frequently occurring FN stage as most representative of the stress state during the period framed by the window, and we associate with this window the a parameter of this FN stage. To overcome the possible influence of the overlapping window approach on the statistical inference, we reduce the range of overlap www.nature.com/scientificreports www.nature.com/scientificreports/ to a maximum of 10 events since the number of nonoverlapping windows is too small to perform the analysis. In this case, the entire series of 509 events is used to construct 12 consecutive windows of 50 elements each with overlaps of 10 events. The statistical test confirms that higher values of the a parameter are associated with lower values of AR, while the FN stages with a <0.6 experience higher AR, with the Z statistics value of the U Mann-Whitney test equal to 2.1, at a significance level p 0 035 = . . Moreover, there is a significant correlation between the a parameter and AR for the tested windows equal to -0.8 (Fig. 7).

Discussion and conclusions
The plot of the mean inferred change in fracture length per unit time versus length implies that the fracturing process in TG is controlled by SFG, with some episodes of instability. The FN growth rate depends on the FN length and the stress parameter. The mean value of parameter a is approximately 0 8 . after correction, and the average value of the stress parameter C is approximately . for cycle 2. Subcritical properties vary with fluid injection loading, and fracture growth occurs with the a parameter ranging from 0 4 . to 1 3 . after correction. At high injection rates, the a parameter is usually high, and at low to mean injection rates, the a parameter varies from low to high values. The values of the fracture growth parameter obtained here are lower than those derived from laboratory tests. Atkinson et al. 2 provided an extensive list of the www.nature.com/scientificreports www.nature.com/scientificreports/ estimated fracture growth index values for Charles' law during laboratory tests. The values of this parameter are found to range from 12 for quartz rocks under water vapor and a temperature of 200 °C to approximately 40 for calcite, basic, granitic and sedimentary rocks and to much higher values for basalts. However, fracture growth index values have also been obtained for geological faults at the crustal scale; they are much lower than those derived from laboratory tests and rarely exceed 4 45 . In our studies, the stress corrosion indices ranges from 0 8 . to . 2 6. These results are similar to those derived for rock mechanical experiments interpreted using the mean field theory for subcritical crack growth based on a modified Griffith criterion for a fractal ensemble of cracks, which range from 0.7 to 3.5 1 . We find that the velocities of fracture growth are much smaller than the seismic rupture propagation. The fracture growth velocity in TG ranges from 10 1 to 10 2 m s , which corresponds to the hydraulic diffusivity in the region of TG on which we focus (e.g., 9 ). This is clear evidence of quasi-static SFG on the system scale, i.e., by intermittent dynamic (seismic) fracture.
Our analysis shows that the proposed scaling for the arrested maximum magnitude with the injected volume of Galis et al. 8 provides good estimates when the injected volume is not high. The temporal evolution of seismicity in TG shows that the observed maximum magnitudes increase the time to the estimated − M w max arr when the injected volume increases. After incorporating the a parameter to estimate the effective bulk modulus for TG in the Galis et al. 8 formula, the maximum observed subcritical magnitudes adequately fit the predicted M w max arr − over the whole time. We relate this behaviour to stress conditions favouring the occurrence of fractures that are randomly distributed in the rock volume and thus decreasing their linkage potential. These conclusions are supported by the results of the analysis of the connectivity process during water injection into the reservoir in TG provided by Orlecka-Sikora et al. 15 and Lasocki and Orlecka-Sikora 46 . The statistical tests confirmed that lower injection rates favoured linking fractures, whereas higher injection rates inhibited linking 15 . This is also in line with the results of our present analysis of the relation between the fracture growth rate and the stress symmetry. Fjaer and Ruisten 40 explained the observed dependency of the rock strength and reduced stress symmetry by linking rock heterogeneity to the impact of the intermediate principal stress. This effect was modelled macroscopically by assigning a smoothened Coulomb failure criterion to each possible orientation of a failure plane in a rock and combining the failure probability of each plane into the overall probability of the failure of the rock. Based on their results, which were obtained at low AR values, only two potential orientations of the failure plane fulfil the Coulomb failure criterion; as a consequence, the expected value of rock strength is higher 40 . We conclude that at lower AR values, the rupturing process may lead to the occurrence of a runaway (critical) rupture because the potential background stress drop of the fault increases. During the analysed period in TG, the strongest seismic event with magnitude Mw 3.2 occurred during the first stage of the second cycle at a mIR value of approximately 6 0 10 m /day 3 3 . ⋅ , where the AR value was low (Fig. 6). The estimated a parameter obtained from the relation presented here at this mIR value is 0.8 (the stress corrosion index exceeds 2).
According to the studies of Fjaer and Ruisten 40 , at high AR values, many orientations of the failure plane are equivalent once the Coulomb failure criterion is fulfilled, and the failure plane assumes the orientation at which the rock fails most easily, e.g., along some rock heterogeneities. This result indicates that fractures may be isolated with low stress drops due to weak rock conditions and that the fracture growth rate decelerates.
Our approach provides a new perspective on earthquake mechanics driven by fluid injection. Recognizing the phenomenon of SFG and its reaction to technological activity and combining this information with the characteristics of the stress state in the reservoir can be used to help manage seismic hazards and optimize technological production. Based on the SFG information, we may conclude that, e.g., n 2 < would imply the transient FN growth rate is slowing 26 , whereas > n 2 would imply increasing event rate and risk of large events during operations. The analysis shows that a higher FN growth rate is related to a lower B value in the Gutenberg-Richter relation (Materials and methods, The relation between b and the fracture growth rate a parameter). This relation can also occur for exponential growth in the expected rate of seismicity with respect to stress or time, as in the Groningen gas field 26,47 , where there have been major interventions in managing the field. The rate of SFG and its relation to the technological stress change factor, which in this case is the water injection rate, provides information about the technological activities that can maintain stable fracture growth, far from a catastrophic runaway earthquake. To better monitor magnitude of seismic events the stress drops and injection information is needed. Seismic hazard management during fluid injection can then be supported by adjusting the injection schedule to keep the fluid pressure and the injected volume below a critical level, at which the observed maximum magnitude approaches the estimated arrested volume at the level set by the local authorities (e.g., 37 ). This insight into the time-and technology-dependent behaviour of rock also helps to identify the stress conditions under which time may be important for the stability of the rock structure, e.g., in mines.

Materials and Methods
Stress drop estimation. Static stress drops were recalculated for the entire input catalogue following the methodology described in 9,10 . Therefore, we provide only a brief description of the source parameter assessment here. The 3-component seismograms obtained from stations located <20 km from the source of seismicity were instrument-corrected and filtered using a 0.2 Hz high-pass Butterworth filter. The 0.5-second time windows with waveforms containing direct P-and S-wave arrivals were tapered using 0.1 s von Hann's taper. The ground P-and S-wave displacement spectra were calculated using the multitaper method 48 separately for each sensor and component. The spectra from each station were then combined 9 . Each spectrum was then fitted separately to the Boatwright point source model, including the frequency-independent attenuation term 9 . The inversion problem, which relies on fitting the curve parameterized by the seismic moment, corner frequency, and quality factor, was optimized by a combination of grid search and simplex techniques. We calculated the median seismic moment and corner frequency for all stations fulfilling the quality criteria. The corresponding source radius was calculated www.nature.com/scientificreports www.nature.com/scientificreports/ assuming the circular source model of Madariaga 49 and a standard rupture velocity of 0.9 V S . Finally, the static stress drop was calculated following Eshelby's formula 14 . Figure 8 presents the range of the earthquake stress drops in the stages of FN growth and a parameter. fracture growth parameter assessment. Table 1 contains the information about the estimates obtained here, together with the standard errors of the growth exponent a and parameter C for 33 of the 66 FN growth periods associated with the three stages of the two considered injection cycles described in the Results section.
Both analysed variables, i.e., the logarithm of the total length of the FN, which is expressed by the cumulative moment magnitude of seismic events, and the logarithm of the mean velocity of FN growth, which is expressed by the cumulative moment magnitude of seismic events divided by the cumulative interevent time, are linearly dependent during the considered stages of water injection, as presented in Table 1. The estimated values of the FN growth parameter a are significant in almost all FN growth stages because their p-values are below the significance level of 0.05. The only two exceptions are the growth parameter a value obtained in the FN marked as 7, which has a p value of 0.054 during the stage proceeding the peak injection in the second cycle of water injection and that in the FN marked as 5, which has a p value of 0.48 during the stage containing the peak injection in the first cycle of injection. The determination coefficient, R 2 , i.e., the proportion of the variance in the dependent variable that can be predicted from the independent variable, is very high, indicating that the regression models built on the total length of the FN explain almost all of the variance in the mean velocity of FN growth in the water injection stages considered. The values of standard error in the regression models are also very low, which indicates the high precision of the model's predictions. In all other stages of water injection cycles, which are not presented in Table 2, we do not observe any linear dependency between the considered variables.
The differences in SFG parameters obtained when assuming the mean velocity of fracture network growth instead of actual velocity in the SfG model. Here, we assume the dependency of fracture growth velocity on fracture length, as described by 2 . For this purpose, we use the observed magnitudes in all identified FNs to calculate the expected velocities of fracture growth for values of the growth parameter a ranging from 0.3 to 1.3. Then, for a given dataset, we calculate the parameter a following the modified equation given by 3 . We then assess how the proposed approach influences the a value and how this impact is related to the number of observations in the dataset, i.e., the number of segments in the FN. Table 2 presents the results of this analysis.
The modified model of FN growth influences the value of the parameter a. We observe that using lower real values of a leads to a smaller underestimation of a. At the lowest considered value of a 0 3 = . , the underestimation is 12%, while at the highest value of = .
a 1 3, the underestimation reaches 37% of the a value. In light of these results, we expect that the a parameter of the FNs in TG, which varies from 0.31 to 0.81, is underestimated and that its real value is higher than that estimated by the FN growth model 2 . The results presented in Table 2 also show the influence of the number of seismic events in the dataset on the estimated results. In Table 2, we present the results obtained using n values corresponding to the number of seismic events in the growth stages of the identified FNs, which are usually less than 30, and those in the entire FN, which are usually less than 100. However, we also studied much higher values of n up to 5000 seismic events in the dataset. Within the range of n values identified in the fracture arrays of the networks in TG, the estimated values of a may vary up to ±0.03. These variations are much lower than the observed variations in the calculated a parameter corresponding to the changes in the mean injection rate. Even at much larger n values, i.e., up to 5000, the variations in the estimated a values are well below 0.1. www.nature.com/scientificreports www.nature.com/scientificreports/ the relation between b and the fracture growth rate a parameter. We consider the exponential magnitude distribution model, which results from the Gutenberg-Richter relation and reads: , where M is the magnitude of events, b is the Gutenberg-Richter constant and M min is the magnitude completeness. The model's parameter is estimated using the maximum likelihood for discrete magnitude values 50,51 : where M is the sample mean of the considered event magnitudes. The analysis is performed in the following steps: 1. Each FN is divided by the intervals of magnitude data with a values: (i) ≤ . a 0 5 and > . a 0 5, (ii) ≤ . a 0 6 and a 0 6 > . , and (iii) ≤ .
a 0 4 and a 0 7 > . . 2. For the selected groups of data, the exponentiality is tested, and the b value is calculated. 3. Then, the difference in b values for different FN growth periods associated with the same FN is tested.
We test the exponentiality of the magnitude distribution using the Anderson-Darling (AD) test 52 (Table 3). Since the AD test is valid for continuous random variables, the magnitudes are randomized within their round-off interval of length 0.01 by the equation proposed by Lasocki and Papadimitriou 53 :  www.nature.com/scientificreports www.nature.com/scientificreports/ where M rand are the randomized magnitudes M, δM is the magnitude round-off interval, u, is a random value from the uniform distribution (0,1), ⋅ F( ) is the cumulative distribution function (CDF) of magnitudes and ⋅ − F ( ) 1 is the corresponding inverse CDF. The magnitude distribution does not follow the exponential distribution when the tested hypothesis of exponentiality is rejected at the significance level of 0.05.
We perform the U Mann-Whitney test to compare the b values between the stages of the FNs with different growth rates. We compare only those pairs of b values wherein the magnitude distributions for both stages are exponential. The results of the analysis are presented in Table 4. The test rejects similarities between the b values of the stages of FNs with the a parameter less than 0.4 and stages with the a parameter greater than 0.7. The FN stages with a higher rate of growth are characterized by lower values of b, so they are more hazardous.
The influence of fracture network growth on the bulk modulus. The presence of cracks in rocks has a significant influence on elastic properties. The effective properties of a rock can be determined as a function of crack development parameters (e.g., [54][55][56]. Generally, the bulk modulus for the cracked rock (K ) is decreased compared to the uncracked rock modulus (K 0 ). Generally, this relation has the form 38 : where B is a constant and ρ is the crack density. For the non-interactive regime, where each crack is unperturbed by other cracks, and the crack-generated strains are summed up, the Eq. (8) is linearized with respect to ρ, and stiffness is linear in crack density: When ρ → 0 two approximations (8) and (9) are compatible. For the non-interactive regime Walsh 57 derived the relation between the bulk modulus K and the bulk modulus of the uncracked material, K 0 , for the isotropic case of randomly oriented penny-shapes cracks: . The crack density depends on the crack interactions, crack orientations, its growth rate, the presence of the initial nuclei and the nucleation rate, the regime in which cracks grow, dilute (freely growing fractures) or dense (stopped fracture), or if they are fluid-filled or not (e.g., 34,35,55,56 Table 2. The results of the estimation of the fracture network growth parameter a using the modified model of fracture network growth given by equation 3 . The input datasets used for the assessment of the modified model performance are the moment magnitudes of the identified fracture networks in The Geysers and the interevent times calculated using the assumed fracture growth parameter a in Charles' law given by 2 ; n is the size of the dataset. www.nature.com/scientificreports www.nature.com/scientificreports/ where n is the number of cracks in the volume V and l is the mean crack length. An exponentially increasing growth rate of crack (Charles' law) produces an exponentially decreasing distribution length (e.g., 35 ): a where a is the growth exponent.  Table 3. Results of the hypothesis for the exponentiality of the magnitude distribution testing using the Anderson-Darling test. The columns are the fracture network code, number of data ( ≤ a a i ), p value from the AD test (a a i ≤ )*, b value ( ≤ a a i ), number of data (a a i > ), p value from AD test ( > a a i )*, b value ( > a a i ), and δb, i.e., difference: b b a a a a i i − ≤ > ; *when p 0 05 < . , exponentiality is rejected at 0.05 significance.
www.nature.com/scientificreports www.nature.com/scientificreports/ The non-interacting approach is generally considered as a good approximation when the crack density is low. For the high crack densities the interactions among cracks become significant and the non-interaction solutions may be inaccurate. To account for crack interactions, numerous improved approximate schemes have been proposed 57 . However, Kachanov 38 reviewing the state of art of the effective elastic properties of cracked solids pointed out that in deriving the effective elastic modulus of cracked rocks, there are still issues needing to be resolved and the accuracy of approximate schemes should be examined since good accuracy derived for specific crack arrays may not apply to other arrays and loading conditions.
For simplicity here we assume the unit volume V, and ( ) . In our case a parameter is below 1 and this function can be approximated by + a 1 . Finally: We found that the most accurate results are when K K 0 in (13) is approximated by − a 1 : When a parameter approaches 1, the A V / is close to 1. B constant is derived empirically from the laboratory tests. In our case B has to be from the interval . (0; 0 5]. The differences between the observed maximum magnitudes and upper limit derived from Galis et al., 8 and for the estimated effective bulk modulus, expressed by RMSE is 0.4 for a (1 ) − coefficient, and 2.2 for the coefficients estimated from the Eqs. (8) and (9). The overestimations of maximum magnitude based on the Eqs. (8) and (9) can be the results of fracture compliance contribution. Here we consider fracture networks thus we expect the coplanar arrangements of fractures. For such fracture organization the dependence of fracture compliance can be stronger than mean radius cubed 38 . Based on these results we approximate the bulk modulus for the fractured rocks in The Geysers by a linear function of the fracture growth rate.

Data availability
The dataset is available on the IS-EPOS platform of the Core Service Anthropogenic Hazards developed in the framework of infrastructure projects IS-EPOS and the European Plate Observing System Programme (https:// tcs.ah-epos.eu). All data needed to evaluate the conclusions in the paper are presented in this paper and in the Materials and methods. Additional data related to this paper may be requested from the authors.  Table 4. The significance value of the U Mann-Whitney test comparing the b values between the stages of the fracture networks with different growth rates as presented in Table 3. Table contains the information about the test statistic, U statistic, and the asymptotic significance (2-tailed) p value.