Analysis of Microseismicity Framing ML > 2.5 Earthquakes at The Geysers Geothermal Field, California

Preparatory mechanisms accompanying or leading to nucleation of larger earthquakes have been observed at both laboratory and field scales, but conditions favoring the occurrence of observable preparatory processes are still largely unknown. In particular, it remains a matter of debate why some earthquakes occur spontaneously without noticeable precursors as opposed to events that are preceded by an extended failure process. In this study, we have generated new high‐resolution seismicity catalogs framing the occurrence of 20 ML > 2.5 earthquakes at The Geysers geothermal field in California. To this end, a seismicity catalog of the 11 days framing each large event was created. We selected 20 sequences sampling different hypocentral depths and hydraulic conditions within the field. Seismic activity and magnitude frequency distributions displayed by the different earthquake sequences are correlated with their location within the reservoir. Sequences located in the northwestern part of the reservoir show overall increased seismic activity and low b values, while the southeastern part is dominated by decreased seismic activity and higher b values. Periods of high injection coincide with high b values and vice versa. These observations potentially reflect varying differential and mean stresses and damage of the reservoir rocks across the field. About 50% of analyzed sequences exhibit no change in seismicity rate in response to the large main event. However, we find complex waveforms at the onset of the main earthquake, suggesting that small ruptures spontaneously grow into or trigger larger events.


Introduction
Characteristic spatial and temporal patterns of seismicity indicating preparation and nucleation of largemagnitude earthquakes would allow to significantly improve short-term earthquake forecasting and subsequent hazard and risk mitigation. For decades, intensive research efforts have been made with the goal of detecting earthquake preparation processes, such as statistical analysis of seismicity catalogs, analysis of foreshocks, or identification of seismic coalescence (e.g., Eneva & Ben-Zion, 1997;Jones, 1985;Mignan, 2012;Mignan, 2014). Since the 1980s, the occurrence of foreshocks accompanied by accelerated seismic moment release (e.g., Jones & Molnar, 1979;Papazachos, 1975;Shaw et al., 1992), increase in the seismicity rates (e.g., Ellsworth et al., 1981), or an increase in Benioff strain (e.g., Ben-Zion & Lyakhovsky, 2002) have been examined as potential signatures to forecast the short-term occurrence of a larger earthquake. Detailed analysis of seismic events involving foreshocks or weak precursor phases as observed for some earthquakes has motivated classical cascade and preslip nucleation models (Ellsworth & Beroza, 1995).
However, the occurrence of such potentially precursory signatures is not systematic in space or time and the conditions favoring them are still a matter of debate (e.g., Kanamori, 1981;Mignan, 2011, and references therein). Some recent large earthquakes have been identified to be preceded by foreshock sequences, such as the 1999 M7.4 Izmit earthquake (Bouchon et al., 2011), the 2010 M7.2 El Mayor-Cucapah earthquake (Hauksson et al., 2010), or the 2011 M9.1 Tohoku-Oki earthquake (e.g., Kato et al., 2012). In contrast, others have shown that foreshocks do not exist systematically (e.g., Zechar & Zhuang, 2010), or that they are indistinguishable from background seismicity (Hardebeck et al., 2008). Furthermore, Jones (1985) estimated that only 6% of the earthquakes in Southern California are followed by an earthquake of larger magnitude within five days and a distance of less than 10 km. This is in contrast to the more recent work of Bouchon et al. (2013) and Brodsky and Lay (2014) who suggest that foreshocks are more common. This is in agreement with many laboratory studies showing foreshock activity prior to stick slip events (Goebel et al., 2012;Lei et al., 2003;Selvadurai & Glaser, 2015). Zaliapin and Ben-Zion (2013a) identified seismicity clusters using nearest-neighbor distributions in a space-time and magnitude domain and found that an increased occurrence of foreshocks occurs in areas displaying high heat flow. In laboratory experiments, acoustic emissions during rock deformation are frequently identified as precursors for the main rupture. A distinct drop in b value preceding the failure of deformed was first described by Scholz (1968) and subsequently confirmed in many studies (Meredith et al., 1990;Zang et al., 1998). Furthermore, Goebel et al. (2012) found that geometric asperities identified in CT scan images were connected to regions of low b values, increased event densities, and moment release over multiple stick-slick cycles. Most recently, Rouet-Leduc et al. (2017) used machine learning to predict the time remaining until failure in shear laboratory experiments by analyzing acoustic emission signals emitted from the fault zone. In fluid-induced seismicity, Walter et al. (2017) identified decreasing b values leading up to the M5.8 Pawnee earthquake in Oklahoma, suggesting that differential stress was slowly building up along the fault in the months prior to the earthquake (Scholz, 2015). Additionally, Savage et al. (2017) observed an increase in localization of foreshocks of the M5.0 Prague earthquake in Oklahoma, initially seen in the entire damage zone and later localizing into a 100-m-thick narrow zone close to the main shock.
In order to fully understand the presence or lack of these precursory signals, we have studied entire earthquake sequences, namely, foreshocks, main shock, aftershocks, and their relation and interplay from The Geysers geothermal field in California. In this context, Martínez-Garzón et al. (2018) previously investigated the relationship between aftershock productivity and injection activity at different geothermal fields in California (including The Geysers) and found a positive correlation of aftershock productivity with netproduced volume. Their results indicate that anthropogenic activity might have a significant influence on foreshock and aftershock patterns in induced seismicity environments.
In addition, high-resolution seismic monitoring allowing for the detection and location of small seismic events may also promote the identification of changing patterns surrounding large event sequences. A recent meta-analysis based on 37 independent studies showed that the nucleation process of earthquakes tends to be detected if the minimum magnitude of the available earthquakes in the catalog is at least 3 orders of magnitude smaller than the selected main shock (Mignan, 2014), typically constraining the detection of earthquake foreshocks (if any) to the largest events included in a seismicity catalog. This is supported by laboratory observations, where high detection capabilities allow anticipating the pending main event before the macroscopic failure occurs (Goebel et al., 2012;Johnson et al., 2013;Kwiatek et al., 2014). Still, several foreshocks signifying the preparation process of a M4.2 earthquake in the eastern Sea of Marmara (Turkey) could be successfully detected using borehole instrumentation (Malin et al., 2018). Yoon et al. (2019) utilized recently developed data mining techniques to investigate the nucleation phase of the 1999 M w 7.1 Hector Mine earthquake, identifying 50 foreshocks in the 20 hr before the main shock. Given already existing data and deployed surface instrumentation, extracting the smallest possible events using state-ofthe-art detection methods is key in identifying physical processes related to large ruptures. That way, Meng et al. (2012) detected~70 times the number of earthquakes listed in the official catalog, when applying a matched filter technique to seismic recordings at the Salton Sea geothermal field surrounding the time of the M w 7.2 El Mayor-Cucapah earthquake. This allowed capturing significant seismicity rate changes due to both dynamic and static stress changes.
The Geysers geothermal field in Northern California shows high seismicity rates, good azimuthal coverage by a local seismic network, and a long history of continuous local high-quality seismic monitoring. The excellent monitoring conditions there allowed collecting one of the most extensive data sets worldwide of fluid-induced seismicity. Here we analyze the microseismicity surrounding the occurrence of 20 moderate-to large-magnitude (2.5 < M L < 4.5) earthquakes at The Geysers geothermal field to investigate how the properties of large earthquake sequences are influenced by ongoing anthropogenic activity and whether large events are preceded by a visible preparation process and how systematic such observations are. To do so, we utilize a matched filter algorithm to extract the maximum amount of data from the continuous seismic recordings from a local high-density broadband network. Merging with existing catalogs and restricting to events located close to the location of each main event, frequency-magnitude distributions are determined aiming to characterize the influence of different geomechanical properties and reservoir characteristics on the behavior of each earthquake sequence. Additionally, a systematic search for potential precursory patterns for each of the sequences is performed. The findings are discussed and related to local stress variations and different damage states within the reservoir rock formation.

Study Region
The Geysers geothermal field is located in the Mayacamas Mountains roughly 110 km north of the San Francisco Bay area in Northern California ( Figure 1). With an installed capacity of more than 1,500 MW across 22 power plants, it is currently the largest geothermal field in the world. Steam is extracted from the vapor-dominated reservoir from approximately 322 active production wells. The field has been producing geothermal energy since the 1960s. Following its peak in 1987, however, production has been declining ever since (Gunasekera et al., 2003), most likely caused by a combination of decreasing reservoir pressure and cooling (Mossop & Segall, 1997). To stabilize production and achieve a net fluid balance within the reservoir, large volumes of treated wastewater are injected across some 54 wells (Majer & Peterson, 2007).
Tectonically, this region is dominated by right-lateral strike-slip motion accommodating the relative movement between the North American Plate and the Pacific Plate (Lisowski et al., 1991). The geothermal field is bounded by two regional faults (Figure 1). Based on GPS-derived slip rates below 1 mm/year the Collayomi Fault is considered to be currently inactive, while the Maacama Fault is considered active based on an average slip rate of about 13 mm/year (Murray et al., 2014). The reservoir is not capped by a geologically well-defined formation, but rather by silica deposition along fluid flow paths resulting in self-sealing of the hydrothermal system (Facca & Tonani, 1965). Today's reservoir is well explored through the large number of wells drilled, and is primarily composed of low porosity (1-2%) but highly fractured greywacke (Lipman et al., 1978). Temperatures across the field vary significantly, between 240°C at 2-km depth in the southeastern and more than 350°C at 2.75 km in the northwestern part of the reservoir (Rutqvist et al., 2016).
During the first half of the twentieth century, The Geysers region was characterized by low seismic activity. Since the active exploitation of the field began, an increasing number of earthquakes have been reported (e.g., Eberhart-Phillips & Oppenheimer, 1984;Oppenheimer, 1986). A dense local short-period network was deployed in 2003, and ever since allows to capture an annual seismicity rate of more than 4,000 earthquakes above magnitude M > 1 (Figure 1). However, seismic activity is not homogeneous across the field, with variations in focal depths (Johnson et al., 2016;Trugman et al., 2016) and frequency magnitude distributions (Convertito et al., 2012;Trugman et al., 2016) between the northwestern and southeastern part of the field, as well as strong seasonal fluctuations in seismicity rate (Johnson et al., 2016;Trugman et al., 2016). Inversion of seismic data for the local maximum horizontal stress shows a general agreement with the NNE-SSW trending regional tectonic structures (Eberhart-Phillips & Oppenheimer, 1984;Martínez-Garzón et al., 2013). The increased seismicity since the start of production has been associated with the ongoing anthropogenic activity through a number of possible physical mechanisms, such as reduced effective stresses due to changes in pore pressure (Majer & Peterson, 2007;Martínez-Garzón et al., 2014). Additionally, temperature contrasts between the injected water and the hot reservoir rock can lead to thermal fracturing close to the injection points (Jeanne et al., 2014;Segall & Fitzgerald, 1998) or cooling-induced geochemical alterations causing microseismic activity (Allis, 1982). Finally, steam production and the resulting reservoir depletion, observed through increased subsidence at the surface (Gunasekera et al., 2003), could further modify the effective stress within the formation.

Data
Seismicity across The Geysers geothermal field is recorded by a permanent surface network of 31 short-period stations operating since 2003, run by the Lawrence Berkeley National Laboratory (Berkeley-Geysers seismic network, IRIS code BG). In the frame of the European Commission GEISER project (EC grant number 241321), this network was supplemented with 26 broadband stations between June 2012 and July 2013, all collocated with selected short-period stations. Each site was equipped with 60-s Guralp or 120-s Trillium sensors, sampled at 200 Hz. The recorded broadband, high-density, and good azimuthal coverage continuous recordings resulted in one of the highest-quality data sets for geothermally induced seismicity. Main noise components were found to be microseism (0.16 Hz) and anthropogenic noise (10-20 Hz; Yu et al., 2018).
Based on automatic processing of the short-period network, two different earthquake catalogs are available through the North California Earthquake Data Center. We utilized the local EGS catalog and corresponding phase picks between June 2012 and July 2013 produced by the Lawrence Berkeley National Laboratory, which represents the most detailed local catalog for the geothermal field (Denllinger et al., 2017). During this time period a total of 120 earthquakes with M L > 3 were observed. Based on this catalog, 20 master events were selected for our analysis ( Figure 2). The locations of these events were chosen as to sample different mean stresses (i.e., by covering the entire possible depth range) as well as fluid pore pressures (i.e., by selecting events occurring during periods of relatively low and high injection flow rates). We selected 20 events starting with the largest magnitudes. If the next potential event was not sampling a new setting of the field, we progressed further down the magnitude list. We chose 10 events located between 2-and 4-km depth, representing the overall reservoir depth range across the majority of the field (although the reservoir tends to be shallower and deeper in the SE and NW, respectively). Out of these 10 events within the reservoir, we took five each during high and low injection periods, respectively. Then, we chose five events from above and five from below the reservoir. During the selection process, we ensured to spatially sample the field as evenly as possible. For each of the 20 master events, between 50 and 120 adjacent events located within their source region, that is, 200-500-m hypocentral distance, were taken as templates for the matched filter technique, supplementing the main shock waveforms. Both P and S wave onsets were manually repicked for each template. Templates and continuous waveforms were detrended and filtered between 4 and 40 Hz using a second-order Butterworth band-pass filter. For each master event and all available stations, we sorted five days of continuous recording before and after the origin day of the master event, together with the day of each master event into 24-hr-long sections.

Matched Filter Algorithm Implementation
The core of the matched filter analysis is a network wide cross correlation between a time series and a known reference signal (e.g., Shearer, 1994;Shelly et al., 2007), where c i (t) is the normalized cross-correlation coefficient at the ith station between the continuous data u i (t) and the template waveform v i (t) (Gibbons & Ringdal, 2006). For big data sets with small sampling intervals this task is computationally expensive, usually requiring a high-performance computing cluster (Beaucé et al., 2017;Meng et al., 2012). Fortunately, the base problem exhibits parallelism on multiple levels, which can significantly reduce computation time (Meng et al., 2012). For each component, multiple templates need to be correlated with the continuous recordings on all stations and components, which can be split up into multiple subseries. Parallelizing all these steps is an obvious way of reducing computational time when memory allocation and sharing is considered.
We implemented the matched filter analysis as follows: (1) a single data input operation is performed that loads the continuous and template data into the memory as separate arrays; (2) the autocorrelation sum of the continuous data and template is precomputed once for later use in the normalization of the correlation calculation (denominator in equation (1)); (3) the correlation is parallelized for each component and template into a maximum number of threads; (4) for each thread the calculation of the cross-correlogram, that is, the correlation coefficient with time, is performed in the frequency domain, as a convolution of the continuous data with the time reversed template; and (5) the cross correlogram is normalized by the precomputed values from (2) ( Figure 3).
After accounting for the station moveout the cross correlograms are stacked across the stations for each template, and positive peaks on the stacked cross correlogram now represent detections with approximately the Figure 2. Location of the 20 master events with M > 2.5 analyzed in this study, as well as the broadband station configuration (green triangles), and the local faults (brown lines). The geothermal field outline is shown by the dashed area, and the dipping reservoir is indicated by the changing shades. Event projections onto the cross section A-A′ are plotted in the bottom left. Color codes are used to distinguish the different earthquake groups: red and blue indicate events located within the reservoir during low and high injection periods, respectively, and yellow and brown indicate events located above and below the reservoir, respectively. same location as the used template event. Detections are extracted by applying a threshold of 7.5 times the median absolute deviation. For a normally distributed random variable the probability of exceeding 7.5 times the median absolute deviation is 2.1 × 10 −7 . Testing the use of different median absolute deviation thresholds and component combinations on small test data sets showed that using only the P wave on vertical component of each station and a template time interval of about 1 s prior to and 3 s after the phase arrival yielded the most robust detection results. We allowed only for a single detection within a 3-s window; in case of multiple detections we kept the detection with the highest mean correlation coefficient across the stations and rejected the remaining.

Phase Picking and Hypocenter Determination
Determination of the exact arrival times for each new detection is based on cross-correlation lag time. We extracted 10-s-long waveforms framing the detection time and cross-correlated them once more with the data from each template at the different stations and components in a narrower time window based on expected travel times along the network, that is, 0.3 s before and 1 s after the manual pick. Vertical components were used to determine the P wave arrival, horizontal components for the S wave. For the template showing the highest cross correlation, the picks were computed as the time lag between the cross correlation maximum and the template pick. We constrained further processing to detections that have at least six picks with at least two picks for each phase (P and S). For each sequence all template phase onsets were first automatically repicked using the correlation with all remaining templates. The new onsets were used to locate each template again. If the location deviated horizontally for more than 150 m from the catalog location that particular template was excluded from the correlation with the detections. Hypocenter locations of events were determined using maximum intersection method, where the earthquake hypocenter is determined independently of the origin time through the use of the equal differential time surface (Font et al., 2004). For any pair of arrivals, such as the same phase on two different stations or two different phases on the same station, an equal differential time surface is defined as the set of all points in the subsurface that satisfy the time difference between the observed arrivals, assuming a fixed velocity model (Zhou, 1994). The intersection of all possible equal differential time surfaces is equal to the global minimum of the cost function used: where δT obs,i is the travel time difference between the observed ith arrival pair, δT th,i is the corresponding theoretical difference, and N is the total number of arrival pairs. Theoretical arrival times were calculated using a local one-dimensional velocity model (Eberhart-Phillips & Oppenheimer, 1984). In the following, we used the Metropolis-Hastings random walk algorithm (Hastings, 1970) to sample the three-dimensional space of hypocenter locations. The standard deviation for the probability density function for each location is on average less than 200 m ( Figure 4). Horizontal location differences for events already located are usually within 150 m, but vertical locations may vary significantly, up to 400 m. It is important to note that the new locations occasionally shift seismic clusters as a whole but do not modify the relative locations of events to each other.

Magnitude of Completeness and b Value Estimation
We calculate the magnitude of each new detected event by comparing the amplitudes between the waveforms and the templates used to find it, assuming that a tenfold increase in amplitude corresponds to 1-unit increase in magnitude: where M temp is the template magnitude and A det and A temp are the peak P wave amplitudes of detection and template, respectively (Peng & Zhao, 2009;Schaff & Richards, 2014;Vuan et al., 2018). However, equation (3) assumes the same location and wave paths for template and detection. To avoid dispersion due to small differences in location we use a site-specific regression: first performing a linear fit between template magnitude and the logarithm of the maximum amplitude for each station, and then calculating a relative magnitude based on that fit. The standard deviation is taken as a measure of error for the relative amplitude determination.
Assuming a Gutenberg-Richter power law distribution of magnitudes, where N is the cumulative number of earthquakes having magnitudes larger than M and a and b are the constants; we calculate the b and a values as a function of the minimum magnitude, M i , using a maximum likelihood estimate (Bender, 1983). Following the goodness-of-fit test of Wiemer and Wyss (2000), we compute a synthetic distribution of magnitudes using the same b, a, and M i values, and calculate the absolute difference, R, of the number of events in each magnitude bin between the observed and synthetic distribution:

10.1029/2019JB017716
Journal of Geophysical Research: Solid Earth where B i and S i are the observed and predicted cumulative number of events in each magnitude bin, respectively. We define the magnitude of completeness, M C , at R = 90%, that is, the point at which a power law can fit 90% of the seismicity or more. For comparison, M C is also calculated using the maximum curvature of the frequency-magnitude distribution (Wiemer & Katsumata, 1999). To assess the stability of the obtained b values we perform a Monte Carlo experiment, where the magnitudes are assigned an error drawn from a normal distribution with a maximum value determined by the magnitude variations in the relative amplitude calculation. This resampling is performed 10 5 times and provides a reliable probability density function for the frequency-magnitude distribution.

Earthquake Catalogs
Applying the matched filter technique to seismic sequences around the 20 selected master events at The Geysers yielded a significant improvement in the detection threshold compared to the existing high-quality local catalogs. We were able to detect a total of more than 63,000 earthquakes in the 11-day periods surrounding 20 master events analyzed (Bentz et al., 2019; Figure 5 and Table 1). Application of the criteria defined for earthquake locations (see section 4) resulted in 47,333 located earthquakes, representing a 75% success ratio between detections and accurate locations. Most of the events that we were unable to locate were part of a small number of sequences, resulting in much higher location ratio (around 85%) for the remaining sequences. This is visible in the congruence between the frequency magnitude distributions from the existing catalogs and the newly derived catalog for the analyzed time periods (Figure 5). Compared to the 32,099 seismic events already existing in the catalog, our new detections almost doubled the number of events. This increased the number of successfully located events by 50% in total. However, this improvement was not uniform across the different sequences, ranging from more than a twofold increase to an actual decrease in located events compared to the local catalog in three cases (Table 1). This decrease in the number of detections was most likely caused by a limited number of templates available for three sequences located near the outer edges of the geothermal field, which resulted in a reduced waveform variability. Comparison of the frequency-magnitude distributions between the original local catalog and a new merged catalog, including all new detections as well as the already located earthquakes included in the Lawrence Berkeley National Laboratory catalog, showed very similar trends ( Figure 5). Both catalogs fulfill a unimodal Gutenberg-Richter relation very well. The distributions tend to diverge for the low magnitudes, demonstrating the catalog improvement achieved.
The increase in detected and located events is much more distinct at close distances to the master event. For this reason, for further analysis we selected only seismic events occurring within a given radius of the master main shock. Different distance thresholds were tested, ranging from 500 m to over 2 km. Assuming standard scaling relations between magnitude, and source radius (Kwiatek et al., 2011) and an average stress drop on the order of 0.1-10 MPa (Kwiatek et al., 2015), the source region for the largest earthquake analyzed was between 1 and 1.5 km. For consistency, we choose to fix the distance threshold across all sequences to a 2-km radius from the main shock. Following this approach, b values with corresponding uncertainties were recalculated for each sequence, restricted to events located within 2 km of the main shock epicenter ( Figure 5 and Table 2). Our Monte Carlo resampling which accounted for the possibility of magnitude errors produced generally stable b values (Table 2 and Figure S1).
In the following, we analyzed potential relations between the statistical properties of events from the 20 catalogs with reservoir characteristics and mechanical (e.g., mean stresses) and hydraulic parameters (e.g., flow rates). Interestingly, we also observed large spatial differences in total seismic activity, that is, the number of events, across the geothermal field. In the northwestern part of The Geysers, the number of events was significantly higher than in the southeastern part (Figure 6d). In depth section the total seismic activity was the highest within the reservoir, with decreased activity above and below the reservoir (Figure 6e). The relationship between seismic activity and main shock magnitude is ambiguous; however, larger-magnitude main shocks tend to also have a higher number of events in their sequence (Figure 6f).

Evolution
When comparing three horizontally collocated sequences with varying depth but in the same area of The Geysers, we observed a larger b value below the reservoir (Figure 7a). At shallow main shock depth the b value was approximately 0.95, while below the reservoir it reached 1.15. The difference exceeded the Note. The number of locations refers to the detections that passed through our strict quality criteria and were located using the maximum intersection method. Catalog improvement is the ratio between the number of locations and the number of events already in the Lawrence Berkeley National Laboratory Enhanced Geothermal System (LBNL EGS) catalog.

10.1029/2019JB017716
Journal of Geophysical Research: Solid Earth   (Figure 7a). At the same time, the total number of events increased with depth for these collocated sequences, from 1,457 at shallow depth to 3,004 at large main shock depth.
In summary, both b value and seismic activity of each sequence were found to be strongly dependent on the epicentral location of the main shock event within the geothermal field, and slightly dependent on the main shock depth. This suggests that seismicity patterns surrounding large-magnitude events at The Geysers are more governed by their location within the field (i.e., the geological and tectonic setting or possibly density of injection wells and production history) than by their depth (i.e., mean stress and temperature).

Changes of b Value and Seismic Activity Related to Injection Volume
The amount of injected water into the reservoir is subjected to strong seasonal fluctuations, allowing us to investigate the influence of injection activity on the behavior of the studied seismic sequences. We identified two sequence pairs, where a large main shock occurred during both high and low injection periods in close proximity to each other (Figures 7b and 7c). At times of high injection volume, an increase in b value could be observed, from 1.15 to 1.20 and from 0.96 to 1.1, respectively. In both cases, the difference in b value is considered significant from the Monte Carlo experiments (Figures 7b and 7c). The largest difference was observed for the sequence pair in the northernmost part of the field, compared to the smaller difference in the central part of the field.
Similarly, the seismic activity varies with the injection volume during the time of the sequence and between different sequence types. Considering the same collocated sequences as discussed above with the b value, we find that for a case located at the northern margin of the investigated area during times of high injection volume about 50% more earthquakes occurred compared to the low injection period. In the central Figure 7. Comparison of b values: (a) between three collocated sequences with different depths and between collocated sequences during high (b) and low (c) injection periods, respectively. In (a), yellow and brown color represent sequences above and below the reservoir, respectively. Blue and red indicate sequences during high and low injection periods, respectively. The location of the sequences on the field is given in the top right corner, as well as the seismicity encoded with hypocenter depths.

10.1029/2019JB017716
Journal of Geophysical Research: Solid Earth portion, this ratio was about 25%. Considering the entire groups of high and low injection sequences, we found a general tendency to increased number of events in each sequence during times when larger fluid volumes were injected ( Figure 6).

Temporal Evolution of b Value Within Each Sequence
We additionally calculated the temporal evolution of b values dividing them into equally sized windows including 100 events. The magnitude of completeness was consistently recalculated for each bin of events. Varying window sizes lead only to smoother or noisier curves but did not affect the general trends. Only seven sequences contained a large enough number of events distributed adequately across time to compute reliable time series of b values, including one above the reservoir. For this latter sequence, a decrease of approximately 0.25 in b value starting three days before the main shock could be observed, with the value stabilizing at a minimum centered on the time of the main event ( Figure 8). Below the reservoir, b values showed a continuous decrease of about 0.2, starting about five days prior to the main shock for each of the two calculated cases. During times of low injection,  (Table 1).

10.1029/2019JB017716
Journal of Geophysical Research: Solid Earth sequences within the reservoir exhibited a relatively stable b value over time, with a small increase of about 0.2, within 24 hr before the main event. In contrast, the two analyzed high injection sequences displayed a similar trend as the shallow sequence, with a decrease of b value starting about three days prior to the main shock ( Figure 8).

Seismicity Rates Surrounding the Main Shock Events
Temporal changes in the seismicity patterns were analyzed by means of the seismicity rate. Equal-sized moving windows containing 10 events were used to calculate the seismicity rate by comparing the origin time of the first and last earthquakes in each bin. Comparing the number of events around the main shocks allowed separating two distinct types of behavior: 11 sequences displayed a significant increase of earthquake rate and seismic moment release after the time of the main event ( Figure 9). However, we did not observe any systematic increase right before the main shock. Almost all sequences within the reservoir during low injection period, and above the reservoir were followed by some sort of seismicity rate increase. In contrast, a second group of events did not show any sign of increased or accelerated seismicity after the time of the main event, those were mostly events during high injection periods or located below the reservoir, that is, times of high stress (Martínez-Garzón et al., 2013). The magnitude of the main shock played an important role, since only events larger than M2.9 showed a significant seismicity increase (Figure 9c). Additionally, the position of each sequence in the field had a visible influence on the rate change pattern. The majority of sequences showing a seismicity rate increase were located in the northwestern part of The Geysers; only a single sequence in the southeastern half displayed a significant rate increase around the time of the main rupture ( Figure 9f). Importantly, changes to seismicity rate were observed only following the main shock, but no sequence displayed an increase of seismicity leading up to the main event. Figure 9. Average earthquake rates around the time of the main shock for each sequence. Earthquake rate (a) 24 hr before and (b) 24 hr following the large-magnitude event compared against main shock magnitude and (d and e) position across the geothermal field (see Figure 2). Yellow diamonds indicate shallow events between 0-and 2-km depth; blue and red show events at 2-4-km depth during high and low injection periods, respectively; and brown represents events below the reservoir. Independent of its magnitude, large-magnitude earthquakes do not always trigger a significant earthquake rate change within the reservoir.

10.1029/2019JB017716
Journal of Geophysical Research: Solid Earth

Discussion
The waveforms from the dense short-period network at The Geysers geothermal field constitute a unique data set of high-quality continuous recordings in the research field of induced seismicity. Extension of the existing permanent short-period network at The Geysers with broadband instrumentation allows the study of earthquakes ranging more than 6 orders of magnitude, from microseismicity to moderately large magnitudes, in full detail. Previous studies on large earthquakes at The Geysers have primarily focused on the relationship between them and anthropogenic activity (e.g., Allis, 1982;Martínez-Garzón et al., 2017;Martínez-Garzón et al., 2018;Ross et al., 1999), its impact on the source physics of earthquakes (e.g., Johnson, 2014;Ross et al., 1996;Yu et al., 2018), and the physical mechanisms governing seismicity (e.g., Kwiatek et al., 2015;Martínez-Garzón et al., 2014). Due to the high resolution and optimal network coverage, seismic data from The Geysers can provide important insights about how reservoir conditions govern the initiation of ruptures.
In the following, we first discuss the resolution of our detection catalogs with respect to other studies. Then, we utilize our observations to discuss whether the large events at The Geysers do display a preparation process and if so, which conditions appear to govern it.

Resolution of the Employed Matched Filter Technique
Application of the matched filter technique to a particular data set may result in a different decrease of the magnitude of completeness depending on the scope of each study. For example, Aiken et al. (2018) investigated potential dynamic triggering at The Geysers by applying a matched filter analysis to continuous vertical recordings half a day before and one day after 25 distant main shocks. The authors were able to lower the magnitude of completeness by 1 full order of magnitude, down to M L~0 , by utilizing around 3,000 templates from the entire geothermal field. The goal in our study was not to increase the overall detection capabilities across the entire geothermal field but to successfully create a detailed record of events occurring in temporal and spatial vicinity of 20 large-magnitude main shocks. This allowed us to capture detailed seismicity patterns, as well as investigate potential preparation processes. Even though the overall magnitude of completeness improvement was only about 0.5 compared to the local catalog, almost 27,000 of our located events constitute previously undetected earthquakes related to the analyzed main shocks. Therefore, direct quantitative comparison between an unbiased local catalog and our detections, which a priori assume a certain waveform, should be done with care.
Our observations highly depend on the quality of the obtained catalogs. The small detection threshold used ensured that we captured events even buried in noise, while at the same time our strict workflow and localization procedure ensured that no regional earthquake was falsely matched and included in our analysis. However, it may have occurred that events were not detected even if they were closely located to the source region of the master events, if their mechanism did not relate to any template mechanism. By selecting a large variety of templates with different waveforms and merging the catalog with the existing local catalog from Lawrence Berkeley National Laboratory this problem should be reduced, but cannot be eliminated entirely.

Sensitivity of Main Shock Sequence Statistics
Seismicity patterns surrounding moderate-magnitude earthquakes at The Geysers are governed by a number of different factors. We found that the position of the main shock within the field plays the largest role in the characteristics displayed by its seismic sequence. The northwestern part of The Geysers shows significantly increased seismic activity around the time of the main shock, compared to the background level. This is reflected in an increase in the daily earthquake rate (Figure 9), as well as an increased total number of events ( Figure 6). On the other hand, large-magnitude main shocks located in the southeastern part generally show no significant change of seismic activity accompanying the main event. Several differences exist between NW and SE of The Geysers, including reservoir depth and the depth of the felsite unit (Jeanne et al., 2014), seismic velocities (Gritto & Jarpe, 2014), coda-Q properties (Blanke et al., 2019), and seismic activity (Johnson et al., 2016), magnitude-frequency statistics of the events (Convertito et al., 2012), and stress field orientation (Martinez-Garzon et al., 2014). In addition, there are more wells injecting fluid in the NW part of the field and the seasonal fluctuations in the volume of water injected are stronger. These features suggest 10.1029/2019JB017716 Journal of Geophysical Research: Solid Earth that reservoir geology and anthropogenic operations might play a dominant role controlling the seismicity associated with moderate-magnitude earthquakes at The Geysers.
The local stress field at The Geysers varies laterally and with depth between normal faulting and strike slip (Martínez-Garzón et al., 2013). However, we did not observe a direct relationship between main shock depth and an accompanying earthquake rate change. The sequences which displayed enhanced seismic activity after the main shock were located across the entire depth range, from 1.7 to 4.2 km, analyzed (Figures 6  and 9), but were restricted to the northwestern part of the field, where seismicity in general and the reservoir itself are located deeper than in the southeast. Additionally, the northwestern part of The Geysers field is characterized by much higher reservoir temperatures, possibly an indication for larger heat flow. Our observation of larger presence of aftershocks in the area of higher temperatures and fluid content is therefore consistent with the findings at global and regional scale from Ben-Zion (2013a, 2016).
Overall, we observe that the main shocks appear to produce increased seismic activity only if the main shock magnitude is M > 2.9, independent of the magnitude of completeness and number of templates. This may indicate that a certain threshold of seismic moment release has to be overcome to transfer sufficient static stress as to trigger aftershocks or to increase seismicity to a significant enough degree to be apparent about the already considerable background rate. Increasing number of triggered aftershocks with main shock magnitude is commonly observed for different seismicity catalogs (Zaliapin & Ben-Zion, 2013b), and it was explained by a combination of the larger magnitude range available to detect aftershocks which are smaller than the main shock and a constant magnitude of completeness. Considering the low magnitude of completeness, however, this cannot explain the complete lack of aftershocks for some of the sequences.
Martínez- Garzón et al. (2013Garzón et al. ( , 2016 showed that the stress field orientation and geometry of reactivated faults additionally vary with time following changes in the fluid injection flow rates. Although seismicity at The Geysers strongly correlate with injection rate and volume, we did not see a strong dependence between the seismic productivity of a main shock and the injected volume around its time of occurrence here ( Figure 6). Similarly, Martínez-Garzón et al. (2018) using cluster analysis found that at The Geysers the number of aftershocks only partly correlates with the injection rates of the field, with the relation varying significantly with time and no continuous trend found.
Experimental results based on rock deformation laboratory experiments have shown that the b value from the frequency magnitude distribution of acoustic events depends inversely on the differential stress (Scholz, 1968(Scholz, , 2015. In agreement with these findings, Schorlemmer et al. (2005) found that b values in Southern California were significantly different for the three different faulting regimes, with larger and lower b values observed for normal and reverse faulting, respectively.
According to poroelasticity, differential stresses should change during time periods of larger fluid injection, particularly along the direction of the maximum horizontal stress (Schoenball et al., 2010), and consequently b values should change. However, we find no average difference between b values observed for sequences during high and low injection volume, respectively. Comparing collocated sequences during different seasons, and thus different injected volumes, we found that during high injection the b value is statistically higher than during times of low injection. These findings are consistent with earlier temporal analysis of b values at one isolated cluster of seismicity in the NW part of the field, where a positive correlation between injection rate and b value over the course of seven years was observed (Leptokaropoulos et al., 2018). Injection-induced increase in pore pressure and the associated reduction of effective stress could lead to weakening of the reservoir rock, resulting in failure before large stresses are able to accumulate. This is in agreement with our observed b value variation across the entire geothermal field, with large b values mostly constrained to the northwestern part, that is, the seismically most active part of The Geysers. Furthermore, the sudden increase of b value at reservoir depth could be explained by the ability to accommodate larger stresses in the unfractured rocks above the geothermal reservoir. Thus, b values at The Geysers might be also influenced by the rock damage within the reservoir (Amitrano, 2003).

Do Preparation Processes Exist for Moderate-Magnitude Earthquakes at The Geysers?
If observed in nature, precursory patterns can be more than 3 orders of magnitude lower than the subsequent main shock (Mignan, 2014), making them extremely difficult to detect and observe for moderate-sized earthquakes. The dense instrumentation at The Geysers provides a low magnitude of completeness data set,

10.1029/2019JB017716
Journal of Geophysical Research: Solid Earth ideal for detecting and studying potential foreshock sequences. Here standard-based energy detectors provide a magnitude of completeness of about M L 1.1, and application of our matched filter reduced further reduced it by half a magnitude for selected sequences, with an even lower detection threshold. However, utilizing these high-quality continuous recordings for 20 large sequences and statistically analyzing microseismicity 4 orders of magnitude lower than the respective main shock, we could not observe any seismic precursory processes in the days preceding the main rupture. None of the sequences analyzed showed an increased spatial correlation of microseismic activity with the main event leading up to its rupture. However, most of the sequences showed localization between subsequent event pairs (see also Figure S1 in the supporting information), but compared with larger groups of events no increased mean correlation could been observed. Decreases of b value leading up to the main failure were seen only for two sequences, thus cannot be generalized. However, similar systematic observations have been made in rock deformation experiments in the laboratory (Goebel et al., 2017;Zang et al., 1998), warranting more detailed analysis of this potential behavior at The Geysers in the future.
Only 11 sequences displayed any reaction in seismic activity to the large-magnitude event, representing an increase of seismicity rate due to aftershock activity and not foreshocks. Similar to these results, Trugman et al. (2016) found anomalously low aftershock productivity when analyzing seismicity across The Geysers, suggesting that earthquakes are primarily induced by localized stress changes from injection, rather than coseismic stresses associated with other earthquakes occurring nearby. In accordance with that, nine of the investigated sequences did not show a response in the seismicity rate to the occurrence of the large-magnitude event. When examining the waveforms for these particular events, smaller events, occurring just seconds before the main failure, were identified on many of them (e.g., Figure 10). Manual picks locate those events within less than a hundred meters of the main shock, indicating that a small area of initial failure propagated to a larger rupture. We observe these complex failures predominantly in the  Table 1).
northeastern and southeastern part of The Geysers. This indicates that in those areas, which have the longest history of injection activity across the field, single small events may propagate into larger ruptures without showing a traditional earthquake sequence with foreshock, main shock, and aftershock. The occurrence of aftershocks was limited almost exclusively to the northwestern part of the geothermal field. Increased injection activity and high temperatures have most likely led to a locally high stress environment, which is much more susceptible to aftershock triggering from stress transfer due to the shorter injection history. Assuming a less damaged formation in this part of the reservoir, the lack of precursory seismic signs of accumulating stress for large ruptures is surprising as the rocks should still be strong enough to display them. This might indicate that stress buildup is not compensated for coseismically but may in fact be aseismic.

Conclusion
We analyzed microseismicity framing the occurrence of 20 moderate-to large-magnitude earthquakes at The Geysers geothermal field, Northern California. Detailed earthquake catalogs for the 11 days surrounding each main shock were obtained using previously detected earthquakes at the nucleation spot as templates in a matched filter algorithm. When considering events located within the rupture plane of each main shock, the magnitude of completeness could be lowered to about M c = 0.5 in most cases. In total approximately 27,000 new unique earthquake locations were obtained.
Applying statistical tools and seismological parameters to the microseismic catalogs, we investigated the seismic behavior before as well as the immediate response of the reservoir following large-magnitude earthquakes. Different parts of the field provide different feedback to the occurrence of large main shocks. In the northwestern part of The Geysers, sequences show increased seismic activity following the main event, while in the southeastern portion of the field no deviation from the background seismicity rate could be observed. Similarly, larger b values for sequences in the northwest than southeast were seen. We interpret this behavior as variations in the local stress field and degree of damage in the reservoir formation. In the northwest high stressing rates result in higher background seismicity rates. Even though no clear relationship between main shock depth and the sequence behavior was observed, the three shallowest sequences exhibit the smallest b values, most likely reflecting the lack of high-damage zones just above the reservoir. Besides the fact of higher density in injection wells in the northwestern The Geysers, no clear temporal link between injection activity and the seismic productivity of each sequence or their b value could be found.
Using statistical analysis, no precursory patterns in the days preceding each large-magnitude event were observed. Seismicity rates prior to main failure were relatively stable, and temporal changes in b value leading up to the main event were seen in only two cases, representing an interesting but isolated observation. Events which did not trigger a change in seismicity rate in the following days show, however, distinct small events in the seconds before their onset, suggesting that in parts of The Geysers large-magnitude events grow from smaller initial earthquakes.
The Geysers reservoir represents a very particular physical environment (e.g., high heat flow, long-lasting fluid injection and production cycles, and shallow reservoir depth). Thus, our results may contribute a new perspective to the ongoing discussion of precursory seismic patterns potentially preceding largemagnitude earthquakes, as well as illuminate specific reservoir responses to stress perturbations dependent on local tectonic and geological settings. Furthermore, we highlight the effectiveness of new earthquake detection and location methodologies on already existing high-quality data sets and the improvement they provide. Additionally, the findings of changing seismic behavior across the geothermal field may have influence on seismic hazard estimations, as they heavily depend on b value und seismicity rate.