Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Evolution of seismicity in relation to fluid injection in the North-Western part of The Geysers geothermal field

Evolution of seismicity in relation to fluid injection in the North-Western part of The Geysers... Summary The Geysers geothermal field located in California, USA, is the largest geothermal site in the world, operating since the 1960s. We here investigate and quantify the correlation between temporal seismicity evolution and variation of the injection data by examination of time-series through specified statistical tools (binomial test to investigate significant rate changes, cross correlation between seismic and injection data, b-value variation analysis). To do so, we utilize seismicity and operational data associated with two injection wells (Prati-9 and Prati-29) which cover a time period of approximately 7 yr (from November 2007 to August 2014). The seismicity is found to be significantly positively correlated with the injection rate. The maximum correlation occurs with a seismic response delay of ∼2 weeks, following injection operations. Those results are very stable even after considering hypocentral uncertainties, by applying a vertical shift of the events foci up to 300 m. Our analysis indicates also time variations of b-value, which exhibits significant positive correlation with injection rates. Hydrothermal systems, Time-series analysis, Induced seismicity, Statistical seismology 1 INTRODUCTION Hundreds of thousands of wells worldwide are used to inject fluid (usually water, occasionally containing chemical additives/contaminants) into the underground for diverse purposes, often causing seismic events occurrence. Well-documented examples of seismic activity induced or triggered by fluid injection are dated back to the 1960s (see Nicholson & Wesson 1990 and Davies et al.2013 for review) and include earthquakes associated with wastewater disposal, secondary oil recovery, solution mining, fluid stimulation to enhance geothermal energy extraction and, during the last decades, hydraulic fracturing for unconventional natural oil and gas exploitation. Such anthropogenic activities may trigger seismic events in critically stressed fault zones, capable to cause failure of pre-existing fractures. In other cases, human operations may entirely control (induce) the nucleation process. Specifically, operations in geothermal fields often lead to considerable increase of seismic activity in areas which were previously characterized by low or even no seismicity and negligible hazard levels (e.g. Majer et al.2007; Mignan et al.2015). The most important feature of seismicity (either induced or triggered) in geothermal fields is its direct correspondence to and dependence on the operational processes, such as fluid injection or production volumes, rates and pressures (e.g. Brodsky & Lajoie 2013). Such relation has been theoretically formulated and established by the non-linear pore pressure diffusion concept (Shapiro & Dinske 2009), or alternatively, by the application of a static stress model based on the non-critical precursory accelerating seismicity theory (Mignan 2016) to induced seismicity. Triggering can generally occur at delay times from hours up to several months since the injection initiation, depending on a variety of factors such as distance from the injection, rock properties, reservoir saturation level, in situ stresses etc. It is however, more usually observed that seismicity onset takes place soon after the initiation of fluid injection (at timescales up to a few days or weeks). For example, Healy et al. (1968) reported a ∼10 d time lag in seismicity following the fluid injection at Rocky Mountain Arsenal. Holland (2013) observed correlations between injection for hydraulic fracturing and seismicity on the scale of a few days, in South-Central Oklahoma. Similarly, cross correlation analysis performed by Kim (2013), between seismicity and injection rate time-series using data from a wastewater disposal well in Youngstown, Ohio, showed that peak seismic activity appears with a time lag of ∼5 d after peak pressures, with the first event detected 13 d after the injection commencement. Evans et al. (2012) reviewed 41 European case histories related to fluid injection in geothermal fields, where in many cases microseismicity starts to be detected some hours after the initiation of operations. The deep understanding of mechanisms involved in the seismogenesis processes should be therefore explored for adjusting injection and production parameters to control seismic activity and mitigate seismic hazard. In order to obtain significant results, high quality data sets should be analysed. For that reason, we selected our study area to be the northwestern part of The Geysers geothermal field, California and utilized a recent relocated seismicity catalog (Martínez-Garzón et al.2014; Kwiatek et al.2015). In the selected area, nearly immediate seismic response to the injection was noticed in several stimulation wells. Based on seismic and operational data from 1976 to 1998 (Majer et al.2007) reported that deep injection-induced seismicity occurs with a short time lag (<2 months). The same authors refer to a particular case where shallow microseismicity was well correlated with injection, rather than production, and occurred with a relatively short time lag of about one week. Stark (2003) performed cross correlation of the monthly seismicity and injection data and obtained a maximum correlation coefficient of 0.52 at a 3 months’ lag providing an interpretation based on rock and water temperature profile. The aim of this study is to investigate the correlation between temporal seismicity evolution (seismic data) and variation of the fluid injection rates (operational data). Correlation analysis is performed and quantified by examination of both seismic and operational data time-series through developed statistical tools. First, we divide the original data set into subsets exhibiting identical either time span or number of events. Then we compare all the subsequent data windows in order to detect and quantify significant variations of seismicity rates and b-values. In doing so, the binomial test is utilized for evaluating the changes in rates of seismic activity and the Aki's (1965) estimator of b-value is used. The derived variations are then compared with the fluctuations of the injection rates in the two nearest injection wells and cross correlation is performed between the time-series of seismic and operational data. In such way, the time delay of seismicity response to the injection process is estimated and the correlation between seismic and operational data is quantified. To strengthen our results, a large range of data windows is tested corresponding to different degrees of data averaging. Hypocentral location (particularly depth) accuracy is also investigated in order to deal with uncertainties rising from the velocity model applied in the relocation process. Our results are combined together in an interpretation mechanism incorporating thermal and pore pressure influence of fluid injection to seismic activity, in agreement with previous numerical and observational models (e.g. Martínez-Garzón et al.2014; Rutqvist et al.2015). 2 STUDY AREA AND DATA The Geysers is the largest producing geothermal field in the world, with many dozens of operating production and injection wells. The thermal energy exploitation takes place since 1960s and the maximum production level was reached in 1987. From that time, the reservoir is recharged by water injections to enhance the steam production and prevent reservoir pressure drawdown. The reservoir is vapour-dominated and built of highly fractured and hydrothermally altered greywacke sandstone (Lockner et al.1982) of low total porosity of about 1–1.5 per cent (Rutqvist et al.2015). It has developed within a complex assemblage of Franciscan rocks ca. 0.25 Ma ago. The temperature of steam is 240 °C at 2 km depth, but exceeds 350 °C in the northwestern part of The Geysers at depths below ca. 2.75 km. No high pressure water injection takes place at The Geysers. Instead, the relatively cool water falls freely into the injection well and into the hot steam reservoir. The steam condenses with a significant volume reduction leading to a negative gauge pressure at the wellhead, which is in contrast to active surface pumping commonly performed for reservoir stimulation with injection at elevated wellhead pressures (e.g. Martínez-Garzón et al.2014). For the analyses performed in the present study, a time-series of 1121 events which occurred up to a maximum distance of 600 m from the open hole of Prati-9 injection well was used (Fig. 1). This seismicity cluster is located in the northwestern part of The Geysers geothermal field and belongs to a data set (Martínez-Garzón et al.2014) consisting of 1254 events which occurred in the vicinity of Prati-9 and Prati-29 injection wells between 2007 December 10 and 2014 of August 23 (Fig. 1). All events were recorded by the Lawrence-Berkeley National Laboratory surface network Berkeley–Geysers (BG) comprising 31 three-component short-period geophones. The events were relocated by Martínez-Garzón et al. (2014) and Kwiatek et al. (2015) from the original Northern California Earthquake Data Center (NCEDC) catalogue, using the double-difference relocation technique (Waldhauser & Ellsworth, 2000). The internal precision of relative hypocentre locations of the events was improved to 50 m due to the applied relocation. Equivalent moment magnitudes were calculated from original duration magnitudes using the formula M*W = 0.9MD + 0.47 (Edwards and Douglas, 2013). The equivalent moment magnitude of completeness as determined by Kwiatek et al. (2015) is M*WC = 1.4. Figure 1. View largeDownload slide Spatial distribution (a—map view, b—cross-section) of all the available data in the study area (1254 events, blue and grey dots). The 1121 seismic events used in our analysis (blue dots) are located in the vicinity of open holes of Prati-9 and Prati-29 injection wells (red circles). Alternative positions of Prati-9 open hole, used for testing the influence of potential events localization vertical shift on the obtained results, are marked with orange circles. Figure 1. View largeDownload slide Spatial distribution (a—map view, b—cross-section) of all the available data in the study area (1254 events, blue and grey dots). The 1121 seismic events used in our analysis (blue dots) are located in the vicinity of open holes of Prati-9 and Prati-29 injection wells (red circles). Alternative positions of Prati-9 open hole, used for testing the influence of potential events localization vertical shift on the obtained results, are marked with orange circles. The selection of this particular data set was done for two reasons. First because of its high quality and small hypocentral uncertainty and second because of the fact that the studied cluster is sufficiently isolated from the neighbouring seismic activity in the broader geothermal field. The selected seismic cluster is concentrated in the vicinity of Prati-9 injection well (average distance of seismic cloud from Prati-9 is 376 m), whereas only a few events (∼30) occurred close to Prati-29 (average distance of seismic cloud from Prati-29 is 800 m). The scarce seismic activity in the close vicinity of Prati-29 led us to exclude those events from our analysis. However, since the fluid volume injected in Prati-29 is proven to influence the total seismic activity (e.g. Kwiatek et al.2015), we considered the effect of the total injection rate and focused on the close vicinity of the open hole of Prati-9 injection well. The injection activity was carried on continuously in Prati-9 (from November 2007 to August 2014) and Prati-29 (from April 2010 to June 2013) wells, with similar average injection rates of 8.79 × 104 m3 month−1 (Kwiatek et al.2015). The injection has seasonal tendency and peak injection rates were observed during the winter months. A pulsation of seismicity cloud was observed following the injection rate changes (Kwiatek et al.2015). Between November 2009 and August 2014, 7.21 and 3.29 Mm3 of treated wastewater were injected into Prati-9 and Prati-29 wells, respectively. 3 METHODOLOGY 3.1 Evaluating the significance of seismicity rate changes—binomial test Let a specified time point be t0. Let also Δt1 and Δt2 be the time intervals and T1, T2 be the time periods before and after t0, defined as [t0 − Δt1, t0] and [t0, t0 + Δt2], respectively. Let the number of events that occurred in the period T1 be n1 and the number of events that occurred in the period T2 be n2. If the seismicity rate during the two periods T1 and T2 is remarkably different, then the actual division of the total number of events in both periods, N = n1 + n2, into n1 and n2 should be significantly different from the division which could be attained at random. Hence the following null hypothesis, H0, is set, stating that n2 could be obtained at random from N under probability P defined as   \begin{equation} P = {{\Delta {t_2}} / {\left( {\Delta {t_2} + \Delta {t_1}} \right)}}. \end{equation} (1) This hypothesis is tested by means of the binomial test (e.g. Wonnacott and Wonnacott, 1977), which provides the probability p1, that if N events occur in a random way in [t0 − Δt1, t0 + Δt2], the number of events in [t0, t0 + Δt2] is less than or equal to n2,   \begin{equation} {p_1} = \Pr \left( {n \le {n_2}|N,P} \right) = \sum\limits_{n = 0}^{{n_2}} {\left( {\begin{array}{@{}*{1}{c}@{}} N\\ n \end{array}} \right){P^n}{{\left( {1 - P} \right)}^{N - n}}} \end{equation} (2)and the probability p2, that the number of events in [t0, t0 + Δt2] is greater than or equal to n2,   \begin{equation} {p_2} = \Pr \left( {n \ge {n_2}\,|\,N,P} \right) = 1 - \sum\limits_{n = 0}^{{n_2} - 1} {\left( {\begin{array}{@{}*{1}{c}@{}} N\\ n \end{array}} \right){P^n}{{\left( {1 - P} \right)}^{N - n}}} . \end{equation} (3) Binomial test treats each event as independent and assigns equal probability of occurrence before and after t0. We test the null hypothesis (H0) at 5 per cent significance level. Therefore if p1 from (2) is less than 0.05, we conclude with the indicated significance that the rate in [t0, t0 + Δt2] decreased with respect to the rate in [t0 – Δt1, t0]. Similarly, if p2 from (3) is less than 0.05, we conclude with the indicated significance that the rate in [t0, t0 + Δ2] increased with respect to the rate in [t0 − Δt1, t0]. 3.2 Binning data Before quantifying the correlation between seismicity and injection we decreased the influence of temporal noise fluctuation on the time-series by averaging (binning) seismic and injection data within specified time windows. We performed binning at different degrees in order to focus on the characteristics of seismicity and injection rates time-series and their correlation in different time scales. Binning was done by averaging the number of events and the values of injection rate over time windows of a selected size. This approach was followed for diminishing the influence of short term, random fluctuations, revealing thus a more robust connection between fluid injection and earthquake occurrence rates. The selection of an optimal time window is crucial at this part of analysis since it may significantly affect the outcome of the process. The window size should be large enough, in order to cover a sufficient time span containing several observations. On the other hand, the window cannot be too large, because it would over-smooth the data. Since there is no objective or commonly accepted criterion determining the optimal degree of binning, we performed trials considering a wide range of window widths, exhibiting either fixed duration of fixed number of data. 3.3 b-value evaluation In order to assess the properties of the size distribution of seismicity we use here the unbounded Gutenberg–Richter (G-R) law, although the actual magnitude distribution of anthropogenic seismicity events may deviate from this distribution (Urban et al.2016). However, as shown in the previously cited work, the deviations are usually small, thus the G-R law can be used as a first approximation of the actual magnitude distribution. The b-value was estimated by the maximum likelihood estimator (Aki, 1965) as follows:   \begin{equation} \hat{b} = \frac{1}{{\ln (10)\left[ {\left\langle M \right\rangle - ({M_{\min }} - \Delta M/2)} \right]}} \end{equation} (4)where 〈M〉 is the sample mean of the events considered, Mmin the completeness threshold and ΔΜ the magnitude binning width. Aki (1965) also provided the estimator of standard deviation of $$\hat{b}$$, σb, defined as:   \begin{equation} {\sigma _b} = \frac{{\hat{b}}}{{\sqrt N }}, \end{equation} (5)where N stands for the sample size. The asymptotic distribution of $$\hat{b}$$ is normal (N($$\hat{b}$$, σb)), and σb can be evaluated by (5), allowing for an adequate estimation of the confidence intervals of $$\hat{b}$$ for relatively large samples. We used this estimator for data sets with more than 100 samples, whereas for smaller data sets bootstrap confidence intervals were calculated instead. To do so, 10 000 realizations of random sampling with replacement were performed for each data set and the confidence intervals were determined corresponding to one standard deviation assuming normal distributed random samples, an assumption that is valid based on the large number of realizations. 4 RESULTS 4.1 Evaluating the significance of seismicity rate changes The results shown in Fig. 2 were obtained by application of the binomial test for identifying considerable seismicity rate changes at 0.05 significance level. For that purpose, we compared time windows of the same duration T, before and after a certain point, t0, which corresponds to T days after the occurrence time of the first event registered in our catalogue. Then we repeated the process by shifting this point by ΔT, and compared the seismicity rate of the corresponding time windows of duration T, before and after this new point, and so on. The time window (T) was selected equal to 50 d and the time step (ΔT) equal to 2 d, thus a total of 1172 data sets were tested. This certain duration of 50 d was selected to ensure that at least one event is included in each one of the windows. The maximum number of earthquakes included in the aforementioned data sets was 58, whereas the maximum difference of events number occurred during subsequent time windows was 40. As shown in Fig. 2, there were significant fluctuations in water injection rates (both in individual injection wells and total) and also in seismicity rates during the entire ∼7 yr's study period. The average rate of ∼0.46 events d−1 exhibits variations, both increases and decreases, from 0.02 to 1.16 events d−1 with the variation of the rate between subsequent time windows being as high as 0.8 events d−1 (for 50 d time windows, overlapping by 2 d). 400 out of the 1172 windows tested were found to exhibit significant seismicity rate changes in comparison with their preceding windows (227 increases and 173 decreases). Similar results were obtained for different window spans of sufficient lengths (10 to 100 d), suggesting that the rate of earthquake occurrence was far from being constant during the period studied. However, it is observed in Fig. 2 that there is a direct association between statistically significant changes of seismicity rates and temporal changes in the flow rate of injected fluid time-series. Fig. 3 shows the comparison between normalized averaged seismicity rates and injection rates for time windows equal to 10, 30 and 50 d. The corresponding scatter plot is demonstrated in Fig. 4 and it verifies the positive correlation between seismicity and operational time-series (Pearson correlation coefficient, rxy, is equal to 0.54, 0.69 and 0.80 for 10, 30 and 50 d windows, respectively). The p-values in all cases are smaller than 10−9 indicating a strong positive correlation between the two averaged time-series which suggests a prompt response of the seismicity rates to changes of injection rate. Some fluctuations are evident (especially for the narrower time window of 10 d), but in general the patterns of these two time-series seem to exhibit considerable mutual similarity, and similar trends. Figure 2. View largeDownload slide Daily total injection rates (blue curve) in relation to the occurrence of earthquakes (stem plot). The daily rates of fluid injected individually into wells Prati-9 and Prati-29 are indicated by black and red curves, respectively. The vertical bars indicate time periods of 50 d which have significantly increased (grey bars) or decreased (green bars) seismicity rates, in comparison with the preceding 50 d window. The time step applied in the analysis is 2 d, so the time windows are overlapping with each other. Figure 2. View largeDownload slide Daily total injection rates (blue curve) in relation to the occurrence of earthquakes (stem plot). The daily rates of fluid injected individually into wells Prati-9 and Prati-29 are indicated by black and red curves, respectively. The vertical bars indicate time periods of 50 d which have significantly increased (grey bars) or decreased (green bars) seismicity rates, in comparison with the preceding 50 d window. The time step applied in the analysis is 2 d, so the time windows are overlapping with each other. Figure 3. View largeDownload slide Comparison between the normalized averaged (binned) time-series: seismicity rates (thick black curve) and total injection rates (blue curve). The dashed grey line denotes the daily total volume of injected fluid. The span of the time windows applied for binning is 10 d (upper frame), 30 d (middle frame) and 50 d (lower frame). Figure 3. View largeDownload slide Comparison between the normalized averaged (binned) time-series: seismicity rates (thick black curve) and total injection rates (blue curve). The dashed grey line denotes the daily total volume of injected fluid. The span of the time windows applied for binning is 10 d (upper frame), 30 d (middle frame) and 50 d (lower frame). Figure 4. View largeDownload slide Scatter plot and Pearson correlation coefficient between binned injection and seismicity rates for different time windows equal to 10 d (left), 30 d (centre) and 50 d (right). Figure 4. View largeDownload slide Scatter plot and Pearson correlation coefficient between binned injection and seismicity rates for different time windows equal to 10 d (left), 30 d (centre) and 50 d (right). 4.2 Correlation between injection and seismicity We calculated the cross correlation function (CCF) between seismicity and injection rates time-series. The cross correlation was performed three times by averaging the injection rates for non-overlapping time windows equal to 1 d (corresponding to the daily measured volume of the injected fluid), 5 and 15 d. The seismicity time-series were obtained as well by calculating the average rate of seismic activity occurred during the aforementioned time windows. For our entire data set of 1121 events the maximum cross correlation value is achieved for a time lag equal to 15 d, meaning that the seismicity follows the injection with a delay of 15 d. In all three approaches followed the correlation was the highest at 15 d time lag, with rxy equal to 0.25, 0.49 and 0.69, for 1 , 5 and 15 d averaged time-series, respectively. The significance of rxy estimate was determined by testing a null hypothesis, according to which, no correlation is evident, that is, rxy essentially equals zero. rxy estimate was compared to 10 000 estimates obtained from randomly selected permutations of the original time-series. If there is no actual causality or correlation between the two original time-series under investigation, then the estimated rxy should not differ significantly from the majority of the values obtained by permutation and therefore the null hypothesis will be true. Otherwise the difference should be large. We selected the significance level for testing the null hypothesis at 0.01 level, meaning than the estimated rxy is higher than at least 99 per cent of the corresponding values arising from the permutated time-series. All of the aforementioned rxy values were found to differ from zero, at 0.01 significance level. Our results indicate that on average there is a short time delay of the order of 2 weeks in the response of seismicity to the fluid injection in the rock matrix. However, as shown in Fig. 5 the time lags corresponding to significant correlation between seismicity and injection rates vary from about −90 to 85 d. If we discard the negative time lags which represent no physical meaning, we obtain an acceptable range of time lags from 0 to ∼85 d, exhibiting a clear peak at 15 d. This suggest that the response of seismic activity to injection operations is more likely to occur with a delay of a few days to ∼3 months, with the strongest correlation calculated at 2-weeks’ time lag. Figure 5. View largeDownload slide CCFs of seismicity and injection data obtained by daily, 5 d and 15 d non-overlapping average (see legend). The x-axis has been limited between −200 and 200 d lag. The rxy values estimated for each time lag are calculated and shown in the y-axis for direct comparison of the graphs. The dashed horizontal lines indicate the 0.99 quantiles obtained in each case from 10 000 data sets, randomly selected with permutation from the corresponding time-series. Figure 5. View largeDownload slide CCFs of seismicity and injection data obtained by daily, 5 d and 15 d non-overlapping average (see legend). The x-axis has been limited between −200 and 200 d lag. The rxy values estimated for each time lag are calculated and shown in the y-axis for direct comparison of the graphs. The dashed horizontal lines indicate the 0.99 quantiles obtained in each case from 10 000 data sets, randomly selected with permutation from the corresponding time-series. The final step was to investigate the influence of shifting the seismicity cloud towards shallower hypocentres. This was done in order to deal with the absolute location uncertainty risen after the relocation process (relative location of the events/velocity model). For this reason, we shifted the injection point 100, 200 and 300 m deeper (equivalent for shifting the events cloud shallower) and recalculated the CCF. Results are shown in Table 1 and they are in good agreement with the ones derived from the original locations. In all depth variation cases studied, rxy are statistically different than zero at 0.01 significance level as well. Table 1. Maximum Pearson correlation coefficient (rxy) for the corresponding cross correlation function value and acceptable time-lags range (1 d averaging of time-series was considered). Results correspond to different absolute hypocentre location of seismicity cloud, shifted along Z axis, from the initial hypocentral depth (Zin). Hypocentre depth  Maximum rxy  Time lag of maximum rxy  Range of accepted time lag values  Initial (Zin)  0.249  15 d  0 to 85 d  Z = Zin + 100 m  0.250  14 d  0 to 91 d  Z = Zin + 200 m  0.252  15 d  0 to 90 d  Z = Zin + 300 m  0.243  19 d  0 to 93 d  Hypocentre depth  Maximum rxy  Time lag of maximum rxy  Range of accepted time lag values  Initial (Zin)  0.249  15 d  0 to 85 d  Z = Zin + 100 m  0.250  14 d  0 to 91 d  Z = Zin + 200 m  0.252  15 d  0 to 90 d  Z = Zin + 300 m  0.243  19 d  0 to 93 d  View Large 4.3 b-value change evaluation In this section, we study the changes of b-value during the ∼7 yr's period of interest. In particular, we investigate the b-value fluctuation in connection with seismicity rate and injection rate changes, according to the following aspects: Distance from the open hole of Prati-9 injection well (spatial b-value distribution). Selected periods considered as stationary in terms of approximately constant seismicity rate (temporal b-value fluctuation for constant seismicity rate). b-value differences in connection with significant seismicity rate differences (temporal b-value fluctuation for overlapping windows exhibiting equal number of events). b-value differences with respect to injection rates. 4.3.1 Distance from the open hole of Prati-9 injection well Table 2 shows the b-values estimated from events which occurred at different hypocentral distances from the open hole section of Prati-9 injection well. The distances were selected in order to form data sets of identical size (∼280 events per set). As shown, the b-value is equal to 1.12 at distances from 300 to 465 m from the open hole of Prati-9, whereas it is slightly higher concerning the near and far field. These differences in b-values are small and they cannot indicate a significant b-value change with the hypocentral distance from the injection open hole. Table 2. b-values and their standard deviations (right column) calculated for different hypocentral distances from Prati-9 (left column). The available sample in each case is shown in the parenthesis. Distance from Prati-9  b-value  $$\phantom{00}$$0–300 m  1.17 ± 0.07 (282 events)  300–380 m  1.12 ± 0.07 (282 events)  380–465 m  1.12 ± 0.07 (278 events)  465–600 m  1.16 ± 0.07 (279 events)  Distance from Prati-9  b-value  $$\phantom{00}$$0–300 m  1.17 ± 0.07 (282 events)  300–380 m  1.12 ± 0.07 (282 events)  380–465 m  1.12 ± 0.07 (278 events)  465–600 m  1.16 ± 0.07 (279 events)  View Large 4.3.2 Selected periods considered as stationary in terms of approximately constant seismicity rate Based on the cumulative occurrence plot (such as Fig. 6), the entire period that the catalog covers was divided visually into 20 approximately stationary parts (straight line segments indicated within consequent different colour-shaded areas, shown in Fig. 6). In each of these parts we estimated the b-value and the standard deviation of its estimate. One can see that b-values significantly vary in time. Some of the data sets include a very small number of events and therefore the corresponding standard deviation bars are large. It is shown that the b-values vary from 0.85 to ∼1.6, a quite large range considering the dimensions of the study area and the time period that the data covers. No correlation of b-value changes with average seismicity rate is found. Figure 6. View largeDownload slide Cumulative seismicity plot (continuous red curve) and 20 segments corresponding to stationary periods (separated by subsequent colour-shaded areas). b-values are calculated for these time periods using eq. (4) (circles) and the corresponding standard deviations are estimated either from eq. (5) (red marks) or as bootstrap 68 per cent confidence intervals (blue marks). The solid horizontal line indicates the average b-value estimated for all available data considering MC = 1.4. Dashed horizontal lines denote ± one standard deviation of the estimated b-value. Total injection rate is represented by solid black line (normalized to its maximum value). Figure 6. View largeDownload slide Cumulative seismicity plot (continuous red curve) and 20 segments corresponding to stationary periods (separated by subsequent colour-shaded areas). b-values are calculated for these time periods using eq. (4) (circles) and the corresponding standard deviations are estimated either from eq. (5) (red marks) or as bootstrap 68 per cent confidence intervals (blue marks). The solid horizontal line indicates the average b-value estimated for all available data considering MC = 1.4. Dashed horizontal lines denote ± one standard deviation of the estimated b-value. Total injection rate is represented by solid black line (normalized to its maximum value). 4.3.3 b-value differences in connection with significant seismicity rate differences In order to study the correlation between changes in earthquake occurrence rate and b-value changes we performed the following analysis. We used a pair of subsequent, non-overlapping windows, each comprising 100 events, that covered the time period in which the first 200 events had been recorded. The constant event window was preferred over the constant time window in order to keep the accuracy of b-value estimation at the same level. We calculated the seismicity rates and b-values in these two windows, respectively, and then we calculated the respective ratios of the results from right-hand-side (later) window to the results from the left-hand-side (former) window. This resulted in two values: the seismicity rate changes ratio and the b-value ratio, which we linked to the time moment of the border between the two windows. We checked the significance of the seismicity rate change using the binomial test. Next we shifted the window pair by 20 events ahead and we repeated the calculations. We continued the procedure until the end of our catalog. The results are demonstrated in Fig. 7, where the Δb-Δrate plot marks the cases of either positive correlation (first and third quadrants) or negative correlation (second and fourth quadrants). No significant correlation between seismicity rate changes and b-value changes can be observed since approximately only 60 per cent of the points are located within the positive correlation regime. We may therefore conclude that seismicity rate changes and b-value variations are uncorrelated. Moreover, there is no obvious systematic temporal pattern in b-value fluctuation. Similar results were obtained for different event windows from 50 to 200 events. No higher than 66 per cent of positive correlation between seismicity rate and b-value changes was achieved in any case. Figure 7. View largeDownload slide Scatter plot with the seismicity rate changes (Δrate) and b-value changes (Δb). Vertical and horizontal dashed lines depict the zero Δrate and zero Δb, respectively. The values were obtained for consequent time windows equal to 50 d, overlapping per 2 d. Positive correlation is evident for approximately 60 per cent of the data sets. The numbers within the figure indicate the quadrants. Figure 7. View largeDownload slide Scatter plot with the seismicity rate changes (Δrate) and b-value changes (Δb). Vertical and horizontal dashed lines depict the zero Δrate and zero Δb, respectively. The values were obtained for consequent time windows equal to 50 d, overlapping per 2 d. Positive correlation is evident for approximately 60 per cent of the data sets. The numbers within the figure indicate the quadrants. 4.3.4 b-value differences with respect to injection rates The b-value analysis so far did not indicate any significant connection between b-value changes and hypocentral distance from Prati-9 injection well open hole (Section 4.3.1), time (Section 4.3.2) and seismic activity rates (Section 4.3.3). Here we directly investigate a potential connection between b-value and injection rate. In doing so, we select time increments of clear either increase or decrease of injection rate and consider the corresponding seismicity. From the several such data sets we chose for our analysis those containing at least 10 events. The results are demonstrated in Fig. 8. With the exception of one data set exhibiting a very high b-value (b = 2.56) all the other data sets have values between 0.73 and 1.60. It is observed that the injection rate changes lead to directly proportional changes of the b-values, with an rxy = 0.57 (excluding the data set with the anomalous b = 2.56). Bootstrap analysis of the injection rate and b-values time-series implies that this correlation is statistically significant at 0.05 level. This means that after randomly resampling the time-series 10 000 times, a higher correlation coefficient was achieved in less than 300 cases (<5 per cent). This result can be also confirmed after stacking the data from decreasing and increasing injection rate periods, respectively, in order to achieve larger and thus, more reliable samples. In this case the b-values correspond to b = 1.12 ± 0.09 (b = 1.18 ± 0.09 including the data set with extreme b-value) and b = 1.33 ± 0.10, for periods of decreasing and increasing injection rate, respectively. Thus, we may conclude that a direct and prompt response between injection rates and b-values is possible. Figure 8. View largeDownload slide b-value estimation for periods of increasing injection rates (red) and decreasing injection rates (blue). The horizontal dashed line indicates the average b-value considering the entire data set. Error bars indicate one standard deviation of the estimated b-values. Figure 8. View largeDownload slide b-value estimation for periods of increasing injection rates (red) and decreasing injection rates (blue). The horizontal dashed line indicates the average b-value considering the entire data set. Error bars indicate one standard deviation of the estimated b-values. 5 DISCUSSION In this section we discuss the overall evolution of seismicity in temporal domain, and review the mechanisms driving the seismicity in this part of the geothermal field. A first observation is that seismic activity is unevenly distributed with distance from the open hole of Prati-9 well (Fig. 9). The characteristic absence of events is shown for the first 100 m from Prati-9 open hole. After the first 100 m, the number of events per unit volume is steadily reduced, with an exception at 300–400 m where an apparent increase of seismic activity is evident. The absence of seismicity in the near field could possibly be linked to the extensive damage caused by the fluid injection resulting in high porosity pathways. The enhancement of seismicity between 300 and 400 m can be explained by the pre-existence of a fault segment which was reactivated by the injection (see Martínez-Garzón et al.2014). Fluid more easily penetrate through such pathways and led to intense seismic activity within a specified zone. Figure 9. View largeDownload slide Number of events per unit volume occurred at different distances from open hole of Prati-9 injection well. Figure 9. View largeDownload slide Number of events per unit volume occurred at different distances from open hole of Prati-9 injection well. A second noteworthy point concerns the propagation of the triggering front and the estimated delay of seismicity. Time delay makes physical sense when interpreted as the time needed for propagation through diffusion of the pressure front to reach particular structures within the reservoir. The pore pressure gradually decreases as distance increases due to the geometrical spreading of the water mass into greater rock volume. According to this, a time lag of even several months is not unexpected and it is controlled by factors such as distance of injection sites from pre-existing fault structures, hydraulic diffusivity and regional stresses (e.g. Brodsky and Lajoie 2013; Langenbruch and Zoback, 2016). In the NW The Geysers area, Johnhson et al. (2016) showed that the seismicity rates during 2005–2007 indicate a 2–5 months’ time lag for all depth intervals, while after 2007 there is no resolvable time lag in the shallow depth intervals (1–2 km), where most of the injection is occurring. Beall et al. (2010) suggested that the delay in elastic response is the time needed to fully saturate portions of the reservoir previously not receiving injection water. Here, in contrast, we study the long term effects of injection (∼7 yr) and not the initial response of the reservoir. Applying typical hydraulic diffusivity values (D = 0.05–0.5 m2 s−1, Talwani et al.2007) to eq. (1) of Shapiro & Dinske (2009), it comes out that seismicity should appear within 600 m distance from injection point within 0.5–6.5 d. If we further consider that the first event was registered 21 d after the injection initiation and was located at ∼500 m from the open hole of Prati-9, this leads to D ≥ 0.011 m2 s−1. The average delay of ∼15 d which we obtained in this analysis can be thus interpreted by the gradual increase of permeability and new fractures generation during the initial phases of injection operations. Such phenomena may accelerate the propagation of both hydraulic and thermal fronts, with the latter considered to advance more slowly (Martínez-Garzón et al.2014; Izadi & Elsworth 2015). Intense seismic activity is concentrated in the close vicinity of Prati-9 (between 100 and 600 m from the open hole), which is located directly within the high temperature zone (HTZ). The main mechanisms driving the seismicity in the analysed area are primarily thermal fracturing and secondary pore pressure diffusion and poroelastic effects (Stark 2003; Martínez-Garzón et al.2014; Rutqvist et al.2015). While thermal stresses are mostly responsible for the majority of seismic events, poroelastic effects may tune the maximum magnitude of seismic activity, especially at larger depths (e.g. Kwiatek et al.2015). The fact that increasing injection rates lead to increased b-values may be interpreted as a result of a sudden increase in pore pressure with the corresponding effective stress reduction (e.g. Scholz 2015). As a result, the crustal strength is decreased, therefore failure occurs before significant stress buildup is accumulated (Staszek et al.2017). Moreover, this leads to a proportional increase of the b-value. This phenomenon has been observed in previous studies (e.g. El-Isa & Eaton 2014 and references therein). It is also worth noting that a strong positive correlation between seismicity rate and injection rate was obtained, whereas no direct correlation between magnitude distribution and seismicity rate was found. This happens because the comparison between b-values and seismicity rates is performed independently of the injection. Seismicity rate changes do not influence b-values neither vice versa, but both parameters are controlled by injection rates. Therefore, disregarding the major cause (i.e. injection) that influences seismic occurrence and magnitude distribution, leads to a failure in revealing any connection between temporal-magnitude seismicity distribution (i.e. only by analysis of overlapping moving time-windows). 6 SUMMARY AND CONCLUSIONS We utilize seismic and operational data from the northwestern part of The Geysers geothermal field (California, USA) to verify and quantify the correlation between injection rates and seismicity rates by analysing the data associated with Prati-9 and Prati-29 injection wells. The correlation between spatio-temporal seismicity evolution and variation of the injection data is performed by examination of original and binned time-series by means of statistical tools (CCF, binomial test for investigating significant rate changes, b-value variation). The major findings of our study are summarized here: A clear positive correlation between seismicity rates and total injection rates is evident. This correlation is statistically different than zero at 0.05 significance level. The entire process seems to take place with a short time lag. We quantified an average time delay of seismic response to injection equal to ∼15 d, whereas a range between 0–85 d is statistically significant at 0.01 level. This time lag is preserved even if the seismicity cloud is shifted at a shallower location as far as 300 m from its original hypocentre. The analysis indicated that b-value changes in time. Moreover, it turned out that b-value is significantly positively correlated with injection rate. The injection rate increases and decreases are followed by increases and decreases of b-value, respectively. No significant influence of the seismicity rates and the distance from the injection well on b-value was detected. ACKNOWLEDGEMENTS This work was supported within SHEER: ‘SHale gas Exploration and Exploitation induced Risks’ project funded from the European Union Horizon 2020 – Research and Innovation Programme, under grant agreement 640896. The work was also partially supported within statutory activities No 3841/E-41/S/2017 of the Ministry of Science and Higher Education of Poland. Original seismic data are available via NCEDC website http://ncedc.org/. The raw injection data and relocated seismic data are available via IS-EPOS platform (https://tcs.ah-epos.eu) after registration and providing affiliation. PMG acknowledges funding from the Helmholtz Association in the frame of their Postdoc Programme. GK acknowledges funding from DFG grant no. KW 84/4–1. REFERENCES Aki K. 1965. Maximum likelihood estimate of b in the formula log N = a − bM and its confidence limits, Bull. Earthq. Res. Inst. Univ. Tokyo , 43, 237– 239. Beall J.J., Wright, M.C., Pingol, A.S., Atkinson P., 2010. Effect of high rate injection on seismicity in The Geysers, Geotherm. Resour. Counc. Trans. , 34, 1203– 1208. Brodsky E.E., Lajoie L.J., 2013. Anthropogenic seismicity rates and operational parameters at the Salton Sea geothermal field, Science , 341( 6145), 543– 546. https://doi.org/10.1126/science.1239213 Google Scholar CrossRef Search ADS PubMed  Davies R., Foulger G., Bindley A., Styles P., 2013. Induced seismicity and hydraulic fracturing for the recovery of hydrocarbons, Mar. Pet. Geol. , 45, 171– 185. https://doi.org/10.1016/j.marpetgeo.2013.03.016 Google Scholar CrossRef Search ADS   Edwards B., Douglas J., 2014. Magnitude scaling of induced earthquakes, Geothermics , 52, 132– 139. https://doi.org/10.1016/j.geothermics.2013.09.012 Google Scholar CrossRef Search ADS   El-Isa Z.H., Eaton D.W., 2014. Spatiotemporal variations in the b-value of earthquake magnitude-frequency distributions: Classification and causes, Tectonophysics , 615–616, 1– 11. https://doi.org/10.1016/j.tecto.2013.12.001 Google Scholar CrossRef Search ADS   Evans K.F., Zappone A., Kraft T., Deichmann N., Moia F., 2012. A survey of the induced seismic responses to fluid injection in geothermal and CO2 reservoirs in Europe, Geothermics , 41, 30– 54. https://doi.org/10.1016/j.geothermics.2011.08.002 Google Scholar CrossRef Search ADS   Healy J.H., Rubey W.W., Griggs D.T., Raleigh C.B., 1968. The Denver earthquakes, Science , 161( 3848), 1301– 1310. https://doi.org/10.1126/science.161.3848.1301 Google Scholar CrossRef Search ADS PubMed  Holland A.A., 2013. Earthquakes triggered by hydraulic fracturing in South-Central Oklahoma, Bull. seism. Soc. Am. , 103( 3), 1784– 1792. https://doi.org/10.1785/0120120109 Google Scholar CrossRef Search ADS   Izadi G., Elsworth D., 2015. The influence of thermal-hydraulic-mechanical- and chemical effects on the evolution of permeability, seismicity and heat production in geothermal reservoirs, Geothermics , 53, 385– 395. https://doi.org/10.1016/j.geothermics.2014.08.005 Google Scholar CrossRef Search ADS   Johnson C.W., Totten E.J., Bürgmann R., 2016. Depth migration of seasonally induced seismicity at The Geysers geothermal field, Geophys. Res. Lett. , 43( 12), 6196– 6204. https://doi.org/10.1002/2016GL069546 Google Scholar CrossRef Search ADS   Kim W., 2013. Induced seismicity associated with fluid injection into a deep well in Youngstown, Ohio, J. geophys. Res. , 118( 7), 3506– 3518. https://doi.org/10.1002/jgrb.50247 Google Scholar CrossRef Search ADS   Kwiatek G., Martínez-Garzón P., Dresen G., Bohnhoff M., Sone H., Hartline C., 2015. Effects of long-term fluid injection on induced seismicity parameters and maximum magnitude in northwestern part of the Geysers geothermal field, J. geophys. Res. , 120( 10), doi:10.1002/2015JB012362. https://doi.org/10.1002/2015JB012362 Langenbruch C., Zoback M.D., 2016. How will induced seismicity in Oklahoma respond to decreased saltwater injection rates?, Sci. Adv. , 2( 11), e1601542, doi:10.1126/sciadv.1601542. https://doi.org/10.1126/sciadv.1601542 Google Scholar CrossRef Search ADS PubMed  Lockner D.A., Summers R., Moore D., Byerlee J.D., 1982. Laboratory measurements of reservoir rock from the Geysers geothermal field, California, Int. J. Rock Mech. Min. Sci. Geomech. Abstr. , 19( 2), 65– 80. https://doi.org/10.1016/0148-9062(82)91632-1 Google Scholar CrossRef Search ADS   Majer E.L., Baria R., Stark M., Oates S., Bommer J., Smith B., Asanuma H., 2007. Induced seismicity associated with enhanced geothermal Systems, Geothermics , 36, 185– 222. https://doi.org/10.1016/j.geothermics.2007.03.003 Google Scholar CrossRef Search ADS   Martínez-Garzón P., Kwiatek G., Sone H., Bohnhoff M., Dresen G., Hartline C., 2014. Spatiotemporal changes, faulting regimes, and source parameters of induced seismicity: A case study from The Geysers geothermal field, J. geophys. Res. , 119, 8378– 8396. https://doi.org/10.1002/2014JB011385 Google Scholar CrossRef Search ADS   Mignan A., 2016. Static behaviour of induced seismicity, Nonlinear Process. Geophys. , 23, 107– 113. https://doi.org/10.5194/npg-23-107-2016 Google Scholar CrossRef Search ADS   Mignan A., Landtwing D., Kästli P., Mena B., Wiemer S., 2015. Induced seismicity risk analysis of the 2006 Basel, Switzerland, Enhanced Geothermal System project: Influence of uncertainties on risk mitigation, Geothermics , 53, 133– 146. https://doi.org/10.1016/j.geothermics.2014.05.007 Google Scholar CrossRef Search ADS   Nicholson C., Wesson R.L., 1990. Earthquake hazard associated with deep well injection; a report to the U.S. Environmental Protection Agency: U.S. Geological Survey Bulletin . Rutqvist J., Dobson P.F., Garcia J., Hartline C., Jeanne P., Oldenburg C.M., Vasco D.W., Walters M., 2015. The northwest geysers EGS demonstration project, California: Pre-stimulation modeling and interpretation of the stimulation, Math. Geosci. , 47, 3– 29. https://doi.org/10.1007/s11004-013-9493-y Google Scholar CrossRef Search ADS   Scholz C.H., 2015. On the stress dependence of the earthquake b value, Geophys. Res. Lett. , 42, doi:10.1002/2014GL062863. https://doi.org/10.1002/2014GL062863 Shapiro S.A., Dinske C., 2009. Scaling of seismicity induced by nonlinear fluid-rock interaction, J. geophys. Res. , 114, 1399, doi: 10.1029/2008JB006145. https://doi.org/10.1029/2008JB006145 Google Scholar CrossRef Search ADS   Stark M., 2003. Seismic evidence for a long-lived enhanced geothermal system (EGS) in the Northern Geysers Reservoir, Geotherm. Resour. Counc. Trans. , 24, 24– 27. Staszek M., Orlecka-Sikora B., Leptokaropoulos K., Kwiatek G., Martínez-Garzón P., 2017. Temporal static stress drop variations due to injection activity at The Geysers geothermal field, California, Geophys. Res. Lett. , 44, doi:10.1002/2017GL073929. https://doi.org/10.1002/2017GL073929 Talwani P., Chen L., Gahalaut K., 2007. Seismogenic permeability, ks, J. geophys. Res. , 112, 7168– 7176. https://doi.org/10.1029/2006JB004665 Google Scholar CrossRef Search ADS   Urban P., Lasocki S., Blascheck P., Do Nascimento A.F., Van Giang N., Kwiatek G., 2016. Violations of Gutenberg–Richter relation in anthropogenic seismicity, Pure appl. Geophys. , doi:10.1007/s00024-015-1188-5. https://doi.org/10.1007/s00024-015-1188-5 Waldhauser F., 2000. A double-difference earthquake location algorithm: method and application to the Northern Hayward Fault, California, Bull. seism. Soc. Am. , 90, 1353– 1368. https://doi.org/10.1785/0120000006 Google Scholar CrossRef Search ADS   Wonnacott T.H., Wonnacott R.J., 1977. Introductory Statistics , Wiley, pp. 650. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

Evolution of seismicity in relation to fluid injection in the North-Western part of The Geysers geothermal field

Download PDF
10 pages

Loading next page...
1
 
/lp/ou_press/evolution-of-seismicity-in-relation-to-fluid-injection-in-the-north-BLuL6if6O4

References (28)

Publisher
Oxford University Press
Copyright
© The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
ISSN
0956-540X
eISSN
1365-246X
DOI
10.1093/gji/ggx481
Publisher site
See Article on Publisher Site

Abstract

Summary The Geysers geothermal field located in California, USA, is the largest geothermal site in the world, operating since the 1960s. We here investigate and quantify the correlation between temporal seismicity evolution and variation of the injection data by examination of time-series through specified statistical tools (binomial test to investigate significant rate changes, cross correlation between seismic and injection data, b-value variation analysis). To do so, we utilize seismicity and operational data associated with two injection wells (Prati-9 and Prati-29) which cover a time period of approximately 7 yr (from November 2007 to August 2014). The seismicity is found to be significantly positively correlated with the injection rate. The maximum correlation occurs with a seismic response delay of ∼2 weeks, following injection operations. Those results are very stable even after considering hypocentral uncertainties, by applying a vertical shift of the events foci up to 300 m. Our analysis indicates also time variations of b-value, which exhibits significant positive correlation with injection rates. Hydrothermal systems, Time-series analysis, Induced seismicity, Statistical seismology 1 INTRODUCTION Hundreds of thousands of wells worldwide are used to inject fluid (usually water, occasionally containing chemical additives/contaminants) into the underground for diverse purposes, often causing seismic events occurrence. Well-documented examples of seismic activity induced or triggered by fluid injection are dated back to the 1960s (see Nicholson & Wesson 1990 and Davies et al.2013 for review) and include earthquakes associated with wastewater disposal, secondary oil recovery, solution mining, fluid stimulation to enhance geothermal energy extraction and, during the last decades, hydraulic fracturing for unconventional natural oil and gas exploitation. Such anthropogenic activities may trigger seismic events in critically stressed fault zones, capable to cause failure of pre-existing fractures. In other cases, human operations may entirely control (induce) the nucleation process. Specifically, operations in geothermal fields often lead to considerable increase of seismic activity in areas which were previously characterized by low or even no seismicity and negligible hazard levels (e.g. Majer et al.2007; Mignan et al.2015). The most important feature of seismicity (either induced or triggered) in geothermal fields is its direct correspondence to and dependence on the operational processes, such as fluid injection or production volumes, rates and pressures (e.g. Brodsky & Lajoie 2013). Such relation has been theoretically formulated and established by the non-linear pore pressure diffusion concept (Shapiro & Dinske 2009), or alternatively, by the application of a static stress model based on the non-critical precursory accelerating seismicity theory (Mignan 2016) to induced seismicity. Triggering can generally occur at delay times from hours up to several months since the injection initiation, depending on a variety of factors such as distance from the injection, rock properties, reservoir saturation level, in situ stresses etc. It is however, more usually observed that seismicity onset takes place soon after the initiation of fluid injection (at timescales up to a few days or weeks). For example, Healy et al. (1968) reported a ∼10 d time lag in seismicity following the fluid injection at Rocky Mountain Arsenal. Holland (2013) observed correlations between injection for hydraulic fracturing and seismicity on the scale of a few days, in South-Central Oklahoma. Similarly, cross correlation analysis performed by Kim (2013), between seismicity and injection rate time-series using data from a wastewater disposal well in Youngstown, Ohio, showed that peak seismic activity appears with a time lag of ∼5 d after peak pressures, with the first event detected 13 d after the injection commencement. Evans et al. (2012) reviewed 41 European case histories related to fluid injection in geothermal fields, where in many cases microseismicity starts to be detected some hours after the initiation of operations. The deep understanding of mechanisms involved in the seismogenesis processes should be therefore explored for adjusting injection and production parameters to control seismic activity and mitigate seismic hazard. In order to obtain significant results, high quality data sets should be analysed. For that reason, we selected our study area to be the northwestern part of The Geysers geothermal field, California and utilized a recent relocated seismicity catalog (Martínez-Garzón et al.2014; Kwiatek et al.2015). In the selected area, nearly immediate seismic response to the injection was noticed in several stimulation wells. Based on seismic and operational data from 1976 to 1998 (Majer et al.2007) reported that deep injection-induced seismicity occurs with a short time lag (<2 months). The same authors refer to a particular case where shallow microseismicity was well correlated with injection, rather than production, and occurred with a relatively short time lag of about one week. Stark (2003) performed cross correlation of the monthly seismicity and injection data and obtained a maximum correlation coefficient of 0.52 at a 3 months’ lag providing an interpretation based on rock and water temperature profile. The aim of this study is to investigate the correlation between temporal seismicity evolution (seismic data) and variation of the fluid injection rates (operational data). Correlation analysis is performed and quantified by examination of both seismic and operational data time-series through developed statistical tools. First, we divide the original data set into subsets exhibiting identical either time span or number of events. Then we compare all the subsequent data windows in order to detect and quantify significant variations of seismicity rates and b-values. In doing so, the binomial test is utilized for evaluating the changes in rates of seismic activity and the Aki's (1965) estimator of b-value is used. The derived variations are then compared with the fluctuations of the injection rates in the two nearest injection wells and cross correlation is performed between the time-series of seismic and operational data. In such way, the time delay of seismicity response to the injection process is estimated and the correlation between seismic and operational data is quantified. To strengthen our results, a large range of data windows is tested corresponding to different degrees of data averaging. Hypocentral location (particularly depth) accuracy is also investigated in order to deal with uncertainties rising from the velocity model applied in the relocation process. Our results are combined together in an interpretation mechanism incorporating thermal and pore pressure influence of fluid injection to seismic activity, in agreement with previous numerical and observational models (e.g. Martínez-Garzón et al.2014; Rutqvist et al.2015). 2 STUDY AREA AND DATA The Geysers is the largest producing geothermal field in the world, with many dozens of operating production and injection wells. The thermal energy exploitation takes place since 1960s and the maximum production level was reached in 1987. From that time, the reservoir is recharged by water injections to enhance the steam production and prevent reservoir pressure drawdown. The reservoir is vapour-dominated and built of highly fractured and hydrothermally altered greywacke sandstone (Lockner et al.1982) of low total porosity of about 1–1.5 per cent (Rutqvist et al.2015). It has developed within a complex assemblage of Franciscan rocks ca. 0.25 Ma ago. The temperature of steam is 240 °C at 2 km depth, but exceeds 350 °C in the northwestern part of The Geysers at depths below ca. 2.75 km. No high pressure water injection takes place at The Geysers. Instead, the relatively cool water falls freely into the injection well and into the hot steam reservoir. The steam condenses with a significant volume reduction leading to a negative gauge pressure at the wellhead, which is in contrast to active surface pumping commonly performed for reservoir stimulation with injection at elevated wellhead pressures (e.g. Martínez-Garzón et al.2014). For the analyses performed in the present study, a time-series of 1121 events which occurred up to a maximum distance of 600 m from the open hole of Prati-9 injection well was used (Fig. 1). This seismicity cluster is located in the northwestern part of The Geysers geothermal field and belongs to a data set (Martínez-Garzón et al.2014) consisting of 1254 events which occurred in the vicinity of Prati-9 and Prati-29 injection wells between 2007 December 10 and 2014 of August 23 (Fig. 1). All events were recorded by the Lawrence-Berkeley National Laboratory surface network Berkeley–Geysers (BG) comprising 31 three-component short-period geophones. The events were relocated by Martínez-Garzón et al. (2014) and Kwiatek et al. (2015) from the original Northern California Earthquake Data Center (NCEDC) catalogue, using the double-difference relocation technique (Waldhauser & Ellsworth, 2000). The internal precision of relative hypocentre locations of the events was improved to 50 m due to the applied relocation. Equivalent moment magnitudes were calculated from original duration magnitudes using the formula M*W = 0.9MD + 0.47 (Edwards and Douglas, 2013). The equivalent moment magnitude of completeness as determined by Kwiatek et al. (2015) is M*WC = 1.4. Figure 1. View largeDownload slide Spatial distribution (a—map view, b—cross-section) of all the available data in the study area (1254 events, blue and grey dots). The 1121 seismic events used in our analysis (blue dots) are located in the vicinity of open holes of Prati-9 and Prati-29 injection wells (red circles). Alternative positions of Prati-9 open hole, used for testing the influence of potential events localization vertical shift on the obtained results, are marked with orange circles. Figure 1. View largeDownload slide Spatial distribution (a—map view, b—cross-section) of all the available data in the study area (1254 events, blue and grey dots). The 1121 seismic events used in our analysis (blue dots) are located in the vicinity of open holes of Prati-9 and Prati-29 injection wells (red circles). Alternative positions of Prati-9 open hole, used for testing the influence of potential events localization vertical shift on the obtained results, are marked with orange circles. The selection of this particular data set was done for two reasons. First because of its high quality and small hypocentral uncertainty and second because of the fact that the studied cluster is sufficiently isolated from the neighbouring seismic activity in the broader geothermal field. The selected seismic cluster is concentrated in the vicinity of Prati-9 injection well (average distance of seismic cloud from Prati-9 is 376 m), whereas only a few events (∼30) occurred close to Prati-29 (average distance of seismic cloud from Prati-29 is 800 m). The scarce seismic activity in the close vicinity of Prati-29 led us to exclude those events from our analysis. However, since the fluid volume injected in Prati-29 is proven to influence the total seismic activity (e.g. Kwiatek et al.2015), we considered the effect of the total injection rate and focused on the close vicinity of the open hole of Prati-9 injection well. The injection activity was carried on continuously in Prati-9 (from November 2007 to August 2014) and Prati-29 (from April 2010 to June 2013) wells, with similar average injection rates of 8.79 × 104 m3 month−1 (Kwiatek et al.2015). The injection has seasonal tendency and peak injection rates were observed during the winter months. A pulsation of seismicity cloud was observed following the injection rate changes (Kwiatek et al.2015). Between November 2009 and August 2014, 7.21 and 3.29 Mm3 of treated wastewater were injected into Prati-9 and Prati-29 wells, respectively. 3 METHODOLOGY 3.1 Evaluating the significance of seismicity rate changes—binomial test Let a specified time point be t0. Let also Δt1 and Δt2 be the time intervals and T1, T2 be the time periods before and after t0, defined as [t0 − Δt1, t0] and [t0, t0 + Δt2], respectively. Let the number of events that occurred in the period T1 be n1 and the number of events that occurred in the period T2 be n2. If the seismicity rate during the two periods T1 and T2 is remarkably different, then the actual division of the total number of events in both periods, N = n1 + n2, into n1 and n2 should be significantly different from the division which could be attained at random. Hence the following null hypothesis, H0, is set, stating that n2 could be obtained at random from N under probability P defined as   \begin{equation} P = {{\Delta {t_2}} / {\left( {\Delta {t_2} + \Delta {t_1}} \right)}}. \end{equation} (1) This hypothesis is tested by means of the binomial test (e.g. Wonnacott and Wonnacott, 1977), which provides the probability p1, that if N events occur in a random way in [t0 − Δt1, t0 + Δt2], the number of events in [t0, t0 + Δt2] is less than or equal to n2,   \begin{equation} {p_1} = \Pr \left( {n \le {n_2}|N,P} \right) = \sum\limits_{n = 0}^{{n_2}} {\left( {\begin{array}{@{}*{1}{c}@{}} N\\ n \end{array}} \right){P^n}{{\left( {1 - P} \right)}^{N - n}}} \end{equation} (2)and the probability p2, that the number of events in [t0, t0 + Δt2] is greater than or equal to n2,   \begin{equation} {p_2} = \Pr \left( {n \ge {n_2}\,|\,N,P} \right) = 1 - \sum\limits_{n = 0}^{{n_2} - 1} {\left( {\begin{array}{@{}*{1}{c}@{}} N\\ n \end{array}} \right){P^n}{{\left( {1 - P} \right)}^{N - n}}} . \end{equation} (3) Binomial test treats each event as independent and assigns equal probability of occurrence before and after t0. We test the null hypothesis (H0) at 5 per cent significance level. Therefore if p1 from (2) is less than 0.05, we conclude with the indicated significance that the rate in [t0, t0 + Δt2] decreased with respect to the rate in [t0 – Δt1, t0]. Similarly, if p2 from (3) is less than 0.05, we conclude with the indicated significance that the rate in [t0, t0 + Δ2] increased with respect to the rate in [t0 − Δt1, t0]. 3.2 Binning data Before quantifying the correlation between seismicity and injection we decreased the influence of temporal noise fluctuation on the time-series by averaging (binning) seismic and injection data within specified time windows. We performed binning at different degrees in order to focus on the characteristics of seismicity and injection rates time-series and their correlation in different time scales. Binning was done by averaging the number of events and the values of injection rate over time windows of a selected size. This approach was followed for diminishing the influence of short term, random fluctuations, revealing thus a more robust connection between fluid injection and earthquake occurrence rates. The selection of an optimal time window is crucial at this part of analysis since it may significantly affect the outcome of the process. The window size should be large enough, in order to cover a sufficient time span containing several observations. On the other hand, the window cannot be too large, because it would over-smooth the data. Since there is no objective or commonly accepted criterion determining the optimal degree of binning, we performed trials considering a wide range of window widths, exhibiting either fixed duration of fixed number of data. 3.3 b-value evaluation In order to assess the properties of the size distribution of seismicity we use here the unbounded Gutenberg–Richter (G-R) law, although the actual magnitude distribution of anthropogenic seismicity events may deviate from this distribution (Urban et al.2016). However, as shown in the previously cited work, the deviations are usually small, thus the G-R law can be used as a first approximation of the actual magnitude distribution. The b-value was estimated by the maximum likelihood estimator (Aki, 1965) as follows:   \begin{equation} \hat{b} = \frac{1}{{\ln (10)\left[ {\left\langle M \right\rangle - ({M_{\min }} - \Delta M/2)} \right]}} \end{equation} (4)where 〈M〉 is the sample mean of the events considered, Mmin the completeness threshold and ΔΜ the magnitude binning width. Aki (1965) also provided the estimator of standard deviation of $$\hat{b}$$, σb, defined as:   \begin{equation} {\sigma _b} = \frac{{\hat{b}}}{{\sqrt N }}, \end{equation} (5)where N stands for the sample size. The asymptotic distribution of $$\hat{b}$$ is normal (N($$\hat{b}$$, σb)), and σb can be evaluated by (5), allowing for an adequate estimation of the confidence intervals of $$\hat{b}$$ for relatively large samples. We used this estimator for data sets with more than 100 samples, whereas for smaller data sets bootstrap confidence intervals were calculated instead. To do so, 10 000 realizations of random sampling with replacement were performed for each data set and the confidence intervals were determined corresponding to one standard deviation assuming normal distributed random samples, an assumption that is valid based on the large number of realizations. 4 RESULTS 4.1 Evaluating the significance of seismicity rate changes The results shown in Fig. 2 were obtained by application of the binomial test for identifying considerable seismicity rate changes at 0.05 significance level. For that purpose, we compared time windows of the same duration T, before and after a certain point, t0, which corresponds to T days after the occurrence time of the first event registered in our catalogue. Then we repeated the process by shifting this point by ΔT, and compared the seismicity rate of the corresponding time windows of duration T, before and after this new point, and so on. The time window (T) was selected equal to 50 d and the time step (ΔT) equal to 2 d, thus a total of 1172 data sets were tested. This certain duration of 50 d was selected to ensure that at least one event is included in each one of the windows. The maximum number of earthquakes included in the aforementioned data sets was 58, whereas the maximum difference of events number occurred during subsequent time windows was 40. As shown in Fig. 2, there were significant fluctuations in water injection rates (both in individual injection wells and total) and also in seismicity rates during the entire ∼7 yr's study period. The average rate of ∼0.46 events d−1 exhibits variations, both increases and decreases, from 0.02 to 1.16 events d−1 with the variation of the rate between subsequent time windows being as high as 0.8 events d−1 (for 50 d time windows, overlapping by 2 d). 400 out of the 1172 windows tested were found to exhibit significant seismicity rate changes in comparison with their preceding windows (227 increases and 173 decreases). Similar results were obtained for different window spans of sufficient lengths (10 to 100 d), suggesting that the rate of earthquake occurrence was far from being constant during the period studied. However, it is observed in Fig. 2 that there is a direct association between statistically significant changes of seismicity rates and temporal changes in the flow rate of injected fluid time-series. Fig. 3 shows the comparison between normalized averaged seismicity rates and injection rates for time windows equal to 10, 30 and 50 d. The corresponding scatter plot is demonstrated in Fig. 4 and it verifies the positive correlation between seismicity and operational time-series (Pearson correlation coefficient, rxy, is equal to 0.54, 0.69 and 0.80 for 10, 30 and 50 d windows, respectively). The p-values in all cases are smaller than 10−9 indicating a strong positive correlation between the two averaged time-series which suggests a prompt response of the seismicity rates to changes of injection rate. Some fluctuations are evident (especially for the narrower time window of 10 d), but in general the patterns of these two time-series seem to exhibit considerable mutual similarity, and similar trends. Figure 2. View largeDownload slide Daily total injection rates (blue curve) in relation to the occurrence of earthquakes (stem plot). The daily rates of fluid injected individually into wells Prati-9 and Prati-29 are indicated by black and red curves, respectively. The vertical bars indicate time periods of 50 d which have significantly increased (grey bars) or decreased (green bars) seismicity rates, in comparison with the preceding 50 d window. The time step applied in the analysis is 2 d, so the time windows are overlapping with each other. Figure 2. View largeDownload slide Daily total injection rates (blue curve) in relation to the occurrence of earthquakes (stem plot). The daily rates of fluid injected individually into wells Prati-9 and Prati-29 are indicated by black and red curves, respectively. The vertical bars indicate time periods of 50 d which have significantly increased (grey bars) or decreased (green bars) seismicity rates, in comparison with the preceding 50 d window. The time step applied in the analysis is 2 d, so the time windows are overlapping with each other. Figure 3. View largeDownload slide Comparison between the normalized averaged (binned) time-series: seismicity rates (thick black curve) and total injection rates (blue curve). The dashed grey line denotes the daily total volume of injected fluid. The span of the time windows applied for binning is 10 d (upper frame), 30 d (middle frame) and 50 d (lower frame). Figure 3. View largeDownload slide Comparison between the normalized averaged (binned) time-series: seismicity rates (thick black curve) and total injection rates (blue curve). The dashed grey line denotes the daily total volume of injected fluid. The span of the time windows applied for binning is 10 d (upper frame), 30 d (middle frame) and 50 d (lower frame). Figure 4. View largeDownload slide Scatter plot and Pearson correlation coefficient between binned injection and seismicity rates for different time windows equal to 10 d (left), 30 d (centre) and 50 d (right). Figure 4. View largeDownload slide Scatter plot and Pearson correlation coefficient between binned injection and seismicity rates for different time windows equal to 10 d (left), 30 d (centre) and 50 d (right). 4.2 Correlation between injection and seismicity We calculated the cross correlation function (CCF) between seismicity and injection rates time-series. The cross correlation was performed three times by averaging the injection rates for non-overlapping time windows equal to 1 d (corresponding to the daily measured volume of the injected fluid), 5 and 15 d. The seismicity time-series were obtained as well by calculating the average rate of seismic activity occurred during the aforementioned time windows. For our entire data set of 1121 events the maximum cross correlation value is achieved for a time lag equal to 15 d, meaning that the seismicity follows the injection with a delay of 15 d. In all three approaches followed the correlation was the highest at 15 d time lag, with rxy equal to 0.25, 0.49 and 0.69, for 1 , 5 and 15 d averaged time-series, respectively. The significance of rxy estimate was determined by testing a null hypothesis, according to which, no correlation is evident, that is, rxy essentially equals zero. rxy estimate was compared to 10 000 estimates obtained from randomly selected permutations of the original time-series. If there is no actual causality or correlation between the two original time-series under investigation, then the estimated rxy should not differ significantly from the majority of the values obtained by permutation and therefore the null hypothesis will be true. Otherwise the difference should be large. We selected the significance level for testing the null hypothesis at 0.01 level, meaning than the estimated rxy is higher than at least 99 per cent of the corresponding values arising from the permutated time-series. All of the aforementioned rxy values were found to differ from zero, at 0.01 significance level. Our results indicate that on average there is a short time delay of the order of 2 weeks in the response of seismicity to the fluid injection in the rock matrix. However, as shown in Fig. 5 the time lags corresponding to significant correlation between seismicity and injection rates vary from about −90 to 85 d. If we discard the negative time lags which represent no physical meaning, we obtain an acceptable range of time lags from 0 to ∼85 d, exhibiting a clear peak at 15 d. This suggest that the response of seismic activity to injection operations is more likely to occur with a delay of a few days to ∼3 months, with the strongest correlation calculated at 2-weeks’ time lag. Figure 5. View largeDownload slide CCFs of seismicity and injection data obtained by daily, 5 d and 15 d non-overlapping average (see legend). The x-axis has been limited between −200 and 200 d lag. The rxy values estimated for each time lag are calculated and shown in the y-axis for direct comparison of the graphs. The dashed horizontal lines indicate the 0.99 quantiles obtained in each case from 10 000 data sets, randomly selected with permutation from the corresponding time-series. Figure 5. View largeDownload slide CCFs of seismicity and injection data obtained by daily, 5 d and 15 d non-overlapping average (see legend). The x-axis has been limited between −200 and 200 d lag. The rxy values estimated for each time lag are calculated and shown in the y-axis for direct comparison of the graphs. The dashed horizontal lines indicate the 0.99 quantiles obtained in each case from 10 000 data sets, randomly selected with permutation from the corresponding time-series. The final step was to investigate the influence of shifting the seismicity cloud towards shallower hypocentres. This was done in order to deal with the absolute location uncertainty risen after the relocation process (relative location of the events/velocity model). For this reason, we shifted the injection point 100, 200 and 300 m deeper (equivalent for shifting the events cloud shallower) and recalculated the CCF. Results are shown in Table 1 and they are in good agreement with the ones derived from the original locations. In all depth variation cases studied, rxy are statistically different than zero at 0.01 significance level as well. Table 1. Maximum Pearson correlation coefficient (rxy) for the corresponding cross correlation function value and acceptable time-lags range (1 d averaging of time-series was considered). Results correspond to different absolute hypocentre location of seismicity cloud, shifted along Z axis, from the initial hypocentral depth (Zin). Hypocentre depth  Maximum rxy  Time lag of maximum rxy  Range of accepted time lag values  Initial (Zin)  0.249  15 d  0 to 85 d  Z = Zin + 100 m  0.250  14 d  0 to 91 d  Z = Zin + 200 m  0.252  15 d  0 to 90 d  Z = Zin + 300 m  0.243  19 d  0 to 93 d  Hypocentre depth  Maximum rxy  Time lag of maximum rxy  Range of accepted time lag values  Initial (Zin)  0.249  15 d  0 to 85 d  Z = Zin + 100 m  0.250  14 d  0 to 91 d  Z = Zin + 200 m  0.252  15 d  0 to 90 d  Z = Zin + 300 m  0.243  19 d  0 to 93 d  View Large 4.3 b-value change evaluation In this section, we study the changes of b-value during the ∼7 yr's period of interest. In particular, we investigate the b-value fluctuation in connection with seismicity rate and injection rate changes, according to the following aspects: Distance from the open hole of Prati-9 injection well (spatial b-value distribution). Selected periods considered as stationary in terms of approximately constant seismicity rate (temporal b-value fluctuation for constant seismicity rate). b-value differences in connection with significant seismicity rate differences (temporal b-value fluctuation for overlapping windows exhibiting equal number of events). b-value differences with respect to injection rates. 4.3.1 Distance from the open hole of Prati-9 injection well Table 2 shows the b-values estimated from events which occurred at different hypocentral distances from the open hole section of Prati-9 injection well. The distances were selected in order to form data sets of identical size (∼280 events per set). As shown, the b-value is equal to 1.12 at distances from 300 to 465 m from the open hole of Prati-9, whereas it is slightly higher concerning the near and far field. These differences in b-values are small and they cannot indicate a significant b-value change with the hypocentral distance from the injection open hole. Table 2. b-values and their standard deviations (right column) calculated for different hypocentral distances from Prati-9 (left column). The available sample in each case is shown in the parenthesis. Distance from Prati-9  b-value  $$\phantom{00}$$0–300 m  1.17 ± 0.07 (282 events)  300–380 m  1.12 ± 0.07 (282 events)  380–465 m  1.12 ± 0.07 (278 events)  465–600 m  1.16 ± 0.07 (279 events)  Distance from Prati-9  b-value  $$\phantom{00}$$0–300 m  1.17 ± 0.07 (282 events)  300–380 m  1.12 ± 0.07 (282 events)  380–465 m  1.12 ± 0.07 (278 events)  465–600 m  1.16 ± 0.07 (279 events)  View Large 4.3.2 Selected periods considered as stationary in terms of approximately constant seismicity rate Based on the cumulative occurrence plot (such as Fig. 6), the entire period that the catalog covers was divided visually into 20 approximately stationary parts (straight line segments indicated within consequent different colour-shaded areas, shown in Fig. 6). In each of these parts we estimated the b-value and the standard deviation of its estimate. One can see that b-values significantly vary in time. Some of the data sets include a very small number of events and therefore the corresponding standard deviation bars are large. It is shown that the b-values vary from 0.85 to ∼1.6, a quite large range considering the dimensions of the study area and the time period that the data covers. No correlation of b-value changes with average seismicity rate is found. Figure 6. View largeDownload slide Cumulative seismicity plot (continuous red curve) and 20 segments corresponding to stationary periods (separated by subsequent colour-shaded areas). b-values are calculated for these time periods using eq. (4) (circles) and the corresponding standard deviations are estimated either from eq. (5) (red marks) or as bootstrap 68 per cent confidence intervals (blue marks). The solid horizontal line indicates the average b-value estimated for all available data considering MC = 1.4. Dashed horizontal lines denote ± one standard deviation of the estimated b-value. Total injection rate is represented by solid black line (normalized to its maximum value). Figure 6. View largeDownload slide Cumulative seismicity plot (continuous red curve) and 20 segments corresponding to stationary periods (separated by subsequent colour-shaded areas). b-values are calculated for these time periods using eq. (4) (circles) and the corresponding standard deviations are estimated either from eq. (5) (red marks) or as bootstrap 68 per cent confidence intervals (blue marks). The solid horizontal line indicates the average b-value estimated for all available data considering MC = 1.4. Dashed horizontal lines denote ± one standard deviation of the estimated b-value. Total injection rate is represented by solid black line (normalized to its maximum value). 4.3.3 b-value differences in connection with significant seismicity rate differences In order to study the correlation between changes in earthquake occurrence rate and b-value changes we performed the following analysis. We used a pair of subsequent, non-overlapping windows, each comprising 100 events, that covered the time period in which the first 200 events had been recorded. The constant event window was preferred over the constant time window in order to keep the accuracy of b-value estimation at the same level. We calculated the seismicity rates and b-values in these two windows, respectively, and then we calculated the respective ratios of the results from right-hand-side (later) window to the results from the left-hand-side (former) window. This resulted in two values: the seismicity rate changes ratio and the b-value ratio, which we linked to the time moment of the border between the two windows. We checked the significance of the seismicity rate change using the binomial test. Next we shifted the window pair by 20 events ahead and we repeated the calculations. We continued the procedure until the end of our catalog. The results are demonstrated in Fig. 7, where the Δb-Δrate plot marks the cases of either positive correlation (first and third quadrants) or negative correlation (second and fourth quadrants). No significant correlation between seismicity rate changes and b-value changes can be observed since approximately only 60 per cent of the points are located within the positive correlation regime. We may therefore conclude that seismicity rate changes and b-value variations are uncorrelated. Moreover, there is no obvious systematic temporal pattern in b-value fluctuation. Similar results were obtained for different event windows from 50 to 200 events. No higher than 66 per cent of positive correlation between seismicity rate and b-value changes was achieved in any case. Figure 7. View largeDownload slide Scatter plot with the seismicity rate changes (Δrate) and b-value changes (Δb). Vertical and horizontal dashed lines depict the zero Δrate and zero Δb, respectively. The values were obtained for consequent time windows equal to 50 d, overlapping per 2 d. Positive correlation is evident for approximately 60 per cent of the data sets. The numbers within the figure indicate the quadrants. Figure 7. View largeDownload slide Scatter plot with the seismicity rate changes (Δrate) and b-value changes (Δb). Vertical and horizontal dashed lines depict the zero Δrate and zero Δb, respectively. The values were obtained for consequent time windows equal to 50 d, overlapping per 2 d. Positive correlation is evident for approximately 60 per cent of the data sets. The numbers within the figure indicate the quadrants. 4.3.4 b-value differences with respect to injection rates The b-value analysis so far did not indicate any significant connection between b-value changes and hypocentral distance from Prati-9 injection well open hole (Section 4.3.1), time (Section 4.3.2) and seismic activity rates (Section 4.3.3). Here we directly investigate a potential connection between b-value and injection rate. In doing so, we select time increments of clear either increase or decrease of injection rate and consider the corresponding seismicity. From the several such data sets we chose for our analysis those containing at least 10 events. The results are demonstrated in Fig. 8. With the exception of one data set exhibiting a very high b-value (b = 2.56) all the other data sets have values between 0.73 and 1.60. It is observed that the injection rate changes lead to directly proportional changes of the b-values, with an rxy = 0.57 (excluding the data set with the anomalous b = 2.56). Bootstrap analysis of the injection rate and b-values time-series implies that this correlation is statistically significant at 0.05 level. This means that after randomly resampling the time-series 10 000 times, a higher correlation coefficient was achieved in less than 300 cases (<5 per cent). This result can be also confirmed after stacking the data from decreasing and increasing injection rate periods, respectively, in order to achieve larger and thus, more reliable samples. In this case the b-values correspond to b = 1.12 ± 0.09 (b = 1.18 ± 0.09 including the data set with extreme b-value) and b = 1.33 ± 0.10, for periods of decreasing and increasing injection rate, respectively. Thus, we may conclude that a direct and prompt response between injection rates and b-values is possible. Figure 8. View largeDownload slide b-value estimation for periods of increasing injection rates (red) and decreasing injection rates (blue). The horizontal dashed line indicates the average b-value considering the entire data set. Error bars indicate one standard deviation of the estimated b-values. Figure 8. View largeDownload slide b-value estimation for periods of increasing injection rates (red) and decreasing injection rates (blue). The horizontal dashed line indicates the average b-value considering the entire data set. Error bars indicate one standard deviation of the estimated b-values. 5 DISCUSSION In this section we discuss the overall evolution of seismicity in temporal domain, and review the mechanisms driving the seismicity in this part of the geothermal field. A first observation is that seismic activity is unevenly distributed with distance from the open hole of Prati-9 well (Fig. 9). The characteristic absence of events is shown for the first 100 m from Prati-9 open hole. After the first 100 m, the number of events per unit volume is steadily reduced, with an exception at 300–400 m where an apparent increase of seismic activity is evident. The absence of seismicity in the near field could possibly be linked to the extensive damage caused by the fluid injection resulting in high porosity pathways. The enhancement of seismicity between 300 and 400 m can be explained by the pre-existence of a fault segment which was reactivated by the injection (see Martínez-Garzón et al.2014). Fluid more easily penetrate through such pathways and led to intense seismic activity within a specified zone. Figure 9. View largeDownload slide Number of events per unit volume occurred at different distances from open hole of Prati-9 injection well. Figure 9. View largeDownload slide Number of events per unit volume occurred at different distances from open hole of Prati-9 injection well. A second noteworthy point concerns the propagation of the triggering front and the estimated delay of seismicity. Time delay makes physical sense when interpreted as the time needed for propagation through diffusion of the pressure front to reach particular structures within the reservoir. The pore pressure gradually decreases as distance increases due to the geometrical spreading of the water mass into greater rock volume. According to this, a time lag of even several months is not unexpected and it is controlled by factors such as distance of injection sites from pre-existing fault structures, hydraulic diffusivity and regional stresses (e.g. Brodsky and Lajoie 2013; Langenbruch and Zoback, 2016). In the NW The Geysers area, Johnhson et al. (2016) showed that the seismicity rates during 2005–2007 indicate a 2–5 months’ time lag for all depth intervals, while after 2007 there is no resolvable time lag in the shallow depth intervals (1–2 km), where most of the injection is occurring. Beall et al. (2010) suggested that the delay in elastic response is the time needed to fully saturate portions of the reservoir previously not receiving injection water. Here, in contrast, we study the long term effects of injection (∼7 yr) and not the initial response of the reservoir. Applying typical hydraulic diffusivity values (D = 0.05–0.5 m2 s−1, Talwani et al.2007) to eq. (1) of Shapiro & Dinske (2009), it comes out that seismicity should appear within 600 m distance from injection point within 0.5–6.5 d. If we further consider that the first event was registered 21 d after the injection initiation and was located at ∼500 m from the open hole of Prati-9, this leads to D ≥ 0.011 m2 s−1. The average delay of ∼15 d which we obtained in this analysis can be thus interpreted by the gradual increase of permeability and new fractures generation during the initial phases of injection operations. Such phenomena may accelerate the propagation of both hydraulic and thermal fronts, with the latter considered to advance more slowly (Martínez-Garzón et al.2014; Izadi & Elsworth 2015). Intense seismic activity is concentrated in the close vicinity of Prati-9 (between 100 and 600 m from the open hole), which is located directly within the high temperature zone (HTZ). The main mechanisms driving the seismicity in the analysed area are primarily thermal fracturing and secondary pore pressure diffusion and poroelastic effects (Stark 2003; Martínez-Garzón et al.2014; Rutqvist et al.2015). While thermal stresses are mostly responsible for the majority of seismic events, poroelastic effects may tune the maximum magnitude of seismic activity, especially at larger depths (e.g. Kwiatek et al.2015). The fact that increasing injection rates lead to increased b-values may be interpreted as a result of a sudden increase in pore pressure with the corresponding effective stress reduction (e.g. Scholz 2015). As a result, the crustal strength is decreased, therefore failure occurs before significant stress buildup is accumulated (Staszek et al.2017). Moreover, this leads to a proportional increase of the b-value. This phenomenon has been observed in previous studies (e.g. El-Isa & Eaton 2014 and references therein). It is also worth noting that a strong positive correlation between seismicity rate and injection rate was obtained, whereas no direct correlation between magnitude distribution and seismicity rate was found. This happens because the comparison between b-values and seismicity rates is performed independently of the injection. Seismicity rate changes do not influence b-values neither vice versa, but both parameters are controlled by injection rates. Therefore, disregarding the major cause (i.e. injection) that influences seismic occurrence and magnitude distribution, leads to a failure in revealing any connection between temporal-magnitude seismicity distribution (i.e. only by analysis of overlapping moving time-windows). 6 SUMMARY AND CONCLUSIONS We utilize seismic and operational data from the northwestern part of The Geysers geothermal field (California, USA) to verify and quantify the correlation between injection rates and seismicity rates by analysing the data associated with Prati-9 and Prati-29 injection wells. The correlation between spatio-temporal seismicity evolution and variation of the injection data is performed by examination of original and binned time-series by means of statistical tools (CCF, binomial test for investigating significant rate changes, b-value variation). The major findings of our study are summarized here: A clear positive correlation between seismicity rates and total injection rates is evident. This correlation is statistically different than zero at 0.05 significance level. The entire process seems to take place with a short time lag. We quantified an average time delay of seismic response to injection equal to ∼15 d, whereas a range between 0–85 d is statistically significant at 0.01 level. This time lag is preserved even if the seismicity cloud is shifted at a shallower location as far as 300 m from its original hypocentre. The analysis indicated that b-value changes in time. Moreover, it turned out that b-value is significantly positively correlated with injection rate. The injection rate increases and decreases are followed by increases and decreases of b-value, respectively. No significant influence of the seismicity rates and the distance from the injection well on b-value was detected. ACKNOWLEDGEMENTS This work was supported within SHEER: ‘SHale gas Exploration and Exploitation induced Risks’ project funded from the European Union Horizon 2020 – Research and Innovation Programme, under grant agreement 640896. The work was also partially supported within statutory activities No 3841/E-41/S/2017 of the Ministry of Science and Higher Education of Poland. Original seismic data are available via NCEDC website http://ncedc.org/. The raw injection data and relocated seismic data are available via IS-EPOS platform (https://tcs.ah-epos.eu) after registration and providing affiliation. PMG acknowledges funding from the Helmholtz Association in the frame of their Postdoc Programme. GK acknowledges funding from DFG grant no. KW 84/4–1. REFERENCES Aki K. 1965. Maximum likelihood estimate of b in the formula log N = a − bM and its confidence limits, Bull. Earthq. Res. Inst. Univ. Tokyo , 43, 237– 239. Beall J.J., Wright, M.C., Pingol, A.S., Atkinson P., 2010. Effect of high rate injection on seismicity in The Geysers, Geotherm. Resour. Counc. Trans. , 34, 1203– 1208. Brodsky E.E., Lajoie L.J., 2013. Anthropogenic seismicity rates and operational parameters at the Salton Sea geothermal field, Science , 341( 6145), 543– 546. https://doi.org/10.1126/science.1239213 Google Scholar CrossRef Search ADS PubMed  Davies R., Foulger G., Bindley A., Styles P., 2013. Induced seismicity and hydraulic fracturing for the recovery of hydrocarbons, Mar. Pet. Geol. , 45, 171– 185. https://doi.org/10.1016/j.marpetgeo.2013.03.016 Google Scholar CrossRef Search ADS   Edwards B., Douglas J., 2014. Magnitude scaling of induced earthquakes, Geothermics , 52, 132– 139. https://doi.org/10.1016/j.geothermics.2013.09.012 Google Scholar CrossRef Search ADS   El-Isa Z.H., Eaton D.W., 2014. Spatiotemporal variations in the b-value of earthquake magnitude-frequency distributions: Classification and causes, Tectonophysics , 615–616, 1– 11. https://doi.org/10.1016/j.tecto.2013.12.001 Google Scholar CrossRef Search ADS   Evans K.F., Zappone A., Kraft T., Deichmann N., Moia F., 2012. A survey of the induced seismic responses to fluid injection in geothermal and CO2 reservoirs in Europe, Geothermics , 41, 30– 54. https://doi.org/10.1016/j.geothermics.2011.08.002 Google Scholar CrossRef Search ADS   Healy J.H., Rubey W.W., Griggs D.T., Raleigh C.B., 1968. The Denver earthquakes, Science , 161( 3848), 1301– 1310. https://doi.org/10.1126/science.161.3848.1301 Google Scholar CrossRef Search ADS PubMed  Holland A.A., 2013. Earthquakes triggered by hydraulic fracturing in South-Central Oklahoma, Bull. seism. Soc. Am. , 103( 3), 1784– 1792. https://doi.org/10.1785/0120120109 Google Scholar CrossRef Search ADS   Izadi G., Elsworth D., 2015. The influence of thermal-hydraulic-mechanical- and chemical effects on the evolution of permeability, seismicity and heat production in geothermal reservoirs, Geothermics , 53, 385– 395. https://doi.org/10.1016/j.geothermics.2014.08.005 Google Scholar CrossRef Search ADS   Johnson C.W., Totten E.J., Bürgmann R., 2016. Depth migration of seasonally induced seismicity at The Geysers geothermal field, Geophys. Res. Lett. , 43( 12), 6196– 6204. https://doi.org/10.1002/2016GL069546 Google Scholar CrossRef Search ADS   Kim W., 2013. Induced seismicity associated with fluid injection into a deep well in Youngstown, Ohio, J. geophys. Res. , 118( 7), 3506– 3518. https://doi.org/10.1002/jgrb.50247 Google Scholar CrossRef Search ADS   Kwiatek G., Martínez-Garzón P., Dresen G., Bohnhoff M., Sone H., Hartline C., 2015. Effects of long-term fluid injection on induced seismicity parameters and maximum magnitude in northwestern part of the Geysers geothermal field, J. geophys. Res. , 120( 10), doi:10.1002/2015JB012362. https://doi.org/10.1002/2015JB012362 Langenbruch C., Zoback M.D., 2016. How will induced seismicity in Oklahoma respond to decreased saltwater injection rates?, Sci. Adv. , 2( 11), e1601542, doi:10.1126/sciadv.1601542. https://doi.org/10.1126/sciadv.1601542 Google Scholar CrossRef Search ADS PubMed  Lockner D.A., Summers R., Moore D., Byerlee J.D., 1982. Laboratory measurements of reservoir rock from the Geysers geothermal field, California, Int. J. Rock Mech. Min. Sci. Geomech. Abstr. , 19( 2), 65– 80. https://doi.org/10.1016/0148-9062(82)91632-1 Google Scholar CrossRef Search ADS   Majer E.L., Baria R., Stark M., Oates S., Bommer J., Smith B., Asanuma H., 2007. Induced seismicity associated with enhanced geothermal Systems, Geothermics , 36, 185– 222. https://doi.org/10.1016/j.geothermics.2007.03.003 Google Scholar CrossRef Search ADS   Martínez-Garzón P., Kwiatek G., Sone H., Bohnhoff M., Dresen G., Hartline C., 2014. Spatiotemporal changes, faulting regimes, and source parameters of induced seismicity: A case study from The Geysers geothermal field, J. geophys. Res. , 119, 8378– 8396. https://doi.org/10.1002/2014JB011385 Google Scholar CrossRef Search ADS   Mignan A., 2016. Static behaviour of induced seismicity, Nonlinear Process. Geophys. , 23, 107– 113. https://doi.org/10.5194/npg-23-107-2016 Google Scholar CrossRef Search ADS   Mignan A., Landtwing D., Kästli P., Mena B., Wiemer S., 2015. Induced seismicity risk analysis of the 2006 Basel, Switzerland, Enhanced Geothermal System project: Influence of uncertainties on risk mitigation, Geothermics , 53, 133– 146. https://doi.org/10.1016/j.geothermics.2014.05.007 Google Scholar CrossRef Search ADS   Nicholson C., Wesson R.L., 1990. Earthquake hazard associated with deep well injection; a report to the U.S. Environmental Protection Agency: U.S. Geological Survey Bulletin . Rutqvist J., Dobson P.F., Garcia J., Hartline C., Jeanne P., Oldenburg C.M., Vasco D.W., Walters M., 2015. The northwest geysers EGS demonstration project, California: Pre-stimulation modeling and interpretation of the stimulation, Math. Geosci. , 47, 3– 29. https://doi.org/10.1007/s11004-013-9493-y Google Scholar CrossRef Search ADS   Scholz C.H., 2015. On the stress dependence of the earthquake b value, Geophys. Res. Lett. , 42, doi:10.1002/2014GL062863. https://doi.org/10.1002/2014GL062863 Shapiro S.A., Dinske C., 2009. Scaling of seismicity induced by nonlinear fluid-rock interaction, J. geophys. Res. , 114, 1399, doi: 10.1029/2008JB006145. https://doi.org/10.1029/2008JB006145 Google Scholar CrossRef Search ADS   Stark M., 2003. Seismic evidence for a long-lived enhanced geothermal system (EGS) in the Northern Geysers Reservoir, Geotherm. Resour. Counc. Trans. , 24, 24– 27. Staszek M., Orlecka-Sikora B., Leptokaropoulos K., Kwiatek G., Martínez-Garzón P., 2017. Temporal static stress drop variations due to injection activity at The Geysers geothermal field, California, Geophys. Res. Lett. , 44, doi:10.1002/2017GL073929. https://doi.org/10.1002/2017GL073929 Talwani P., Chen L., Gahalaut K., 2007. Seismogenic permeability, ks, J. geophys. Res. , 112, 7168– 7176. https://doi.org/10.1029/2006JB004665 Google Scholar CrossRef Search ADS   Urban P., Lasocki S., Blascheck P., Do Nascimento A.F., Van Giang N., Kwiatek G., 2016. Violations of Gutenberg–Richter relation in anthropogenic seismicity, Pure appl. Geophys. , doi:10.1007/s00024-015-1188-5. https://doi.org/10.1007/s00024-015-1188-5 Waldhauser F., 2000. A double-difference earthquake location algorithm: method and application to the Northern Hayward Fault, California, Bull. seism. Soc. Am. , 90, 1353– 1368. https://doi.org/10.1785/0120000006 Google Scholar CrossRef Search ADS   Wonnacott T.H., Wonnacott R.J., 1977. Introductory Statistics , Wiley, pp. 650. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

Journal

Geophysical Journal International – Oxford University Press

Published: Feb 1, 2018

There are no references for this article.