Matched filter earthquake detection relies on waveform cross correlation of a known earthquake signal or signals with continuously recorded seismic data. It has been shown to be especially useful in noisy areas where more traditional, amplitude-based techniques fail \citep{Gibbons_2006, Shelly_2007}. This makes it particularly useful in areas of geothermal power generation where large numbers of highly similar, low-magnitude earthquakes can occur and which are usually subject to high levels of anthropogenic noise associated with the power generation process. We have already completed Objective 1 for the 2015 earthquake catalog and outline the methodology below. The Mercury datasets for 2012-2014 will be processed in the same way.
The 637 events taken from the filtered GNS catalog outlined in Section \ref{Dataset} were used as ’template events’ (Figure \ref{Figure1}). Each template waveform was cut to a length of one second, starting 0.1 seconds before the P-pick on each channel for which a pick was made. Waveforms were filtered from 1.0 to 20.0 Hz after confirming that most events in the dataset contained significant amounts of energy within that pass band. Continuous seismic data were processed in an identical manner to the templates within the matched filter routine.
To generate detections, event templates were cross correlated with continuous data at a rate of 50 samples per second. At each sample, the cross correlation coefficients for each channel of data were summed to create the network detection statistic. Whenever the detection statistic exceeded a threshold value, which in this case was defined as the median absolute deviation (MAD) of the detection statistic multiplied by 8 (as suggested by \citep{Shelly_2007}), a detection was recorded.
After removing duplicate detections, P-picks for each of the newly detected events were made based upon cross correlation of the detecting template with the continuous data at the time of the detection. Template waveforms were correlated with network data over a 0.2 second window centered on the detection time. Picks were recorded at the time corresponding to the highest correlation value within that window. Time errors were then assigned to each pick based upon its correlation value with the template. Through visual review of a random sample of 100 events we assigned timing uncertainties of 0.01 seconds to picks with correlation values of greater than 0.90, 0.05 s for values of between 0.70 and 0.90, 0.10 seconds for values between 0.50 and 0.70, 0.30 s for values between 0.30 and 0.50 and 0.50 seconds for all other picks. We kept all events with more than 7 picks and discarded those with fewer. From an original total of nearly 30,000 detections, we were left with 15,471 detections after removing duplicates and discarding events with an insufficient number of picks.
For these 15,471 events, we then incorporated S-picks made by MSc student Stefan Mroczek \citep{mroczek2016} using the method developed by Diehl and others \citep{Diehl_2009}. These events were then located with the nonlinear location program NonLinLoc \citep{Lomax_2014} using a 1-D model for Rotokawa used by former VUW MSc student Zara Rawlinson (\citep{rawlinson2011microseismicity} pp. 77).
In the future, we will investigate the effect that S-picks have on the location uncertainties at Rotokawa and Ngatamariki. In the past it has been reported that S-picks, when added to a catalog, actually have the effect of increasing location uncertainties (Steve Sherburn, personal communication). This increase in uncertainty is probably a reflection of the extremely heterogeneous local velocity structure at the geothermal fields and that we are grossly oversimplifying this structure by assuming a 1-dimensional velocity model. To explore this further, we will make use of a 3-dimensional model for both fields being developed by PhD student Steve Sewell. Comparing locations using only P-picks with locations using both P and S-picks for both a 1-D and 3-D model should help determine the sensitivity of locations to velocity model variations, particularly the S-wave velocity structure.
As a final step, the entire catalog was relocated using the double-difference relocation program HypoDD \citep{Waldhauser_2000} with differential pick times generated using a cross correlation routine within the EQcorrscan Python package \citep{Chamberlain_2014}. Due to the size of the catalog and memory limitations in the software, events were split into separate Rotokawa and Ngatamariki catalogs before being relocated with HypoDD.
Specific periods from 2012 until the end of 2015 have been outlined by Mercury as having particular significance from a resource management standpoint. These center around large scale changes in injection, for example drops in injectivity at a particular well or shifts in injection to a new well, well stimulation activities (at Ngatamariki, specifically) and drilling of new wells. The periods outlined are the following:
Ngatamariki
NM9 stimulation (December, 2012-January, 2013)
Plant start-up and associated stimulation of injection wells (April-August, 2013)
Total losses during drilling of NM12 (June-September, 2014)
Rotokawa
RK24 injectivity decline (2012-2013)
Any period of ’swarm’ behaviour, specifically during plant startup and shutdown (>10 events per day)
Total losses during drilling of RK34 (September-November, 2014)
Switch of injection to RK22 (August-September, 2015)
To compute the magnitudes of the events detected using the matched filter method outlined above, we will use the SVD method detailed by Rubinstein and Ellsworth \citep{Rubinstein_2010}. This method makes use of the waveform similarity inherent between events in repeating sequences, such as the sequences defined by matched filter detections for a single template event, to calculate the relative moment of all events in the sequence. To implement this, we will separate the matched filter detections out into sequences based first on detecting template, and then on waveform similarity and hypocentral location. As each template event was detected and processed automatically, a local magnitude has already been assigned to each event. We can use these pre-assigned magnitudes to calibrate the relative moment calculations from the SVD method, provided each sequence contains only one template event. In this way we will be left with accurate local magnitude calculations based on relative moment calculations via the SVD method.
Calculation of the focal mechanisms for the matched filter-detected events in this study will be done using the approach outlined in \citep{Pugh_2016} which expresses the first motion of a given phase arrival as the probability of both negative and positive polarities instead of the standard, binary “up or down” approach. This approach can also be fully automated which is useful when dealing with datasets of the size of matched filter detection catalogs. \citep{Pugh_2016, Pugh_2016a}
Match filter detection techniques are computationally expensive when applied across large datasets. The expense is multiplied as more earthquakes are added to the list of desired template events (referred to by Barrett and Beroza as the design set) \citep{Barrett_2014} and as larger seismograph networks, with more channels of data, are used. At the same time, it is important that a design set effectively represent the variety of sources present in a given study area or risk missing what might be important seismicity. Subspace detection has been used to address these considerations (e.g. \citep{Harris_2006, Harris_2006a, Gibbons_2006, Barrett_2014}). Subspace detection tries to represent a design set of template events, which would be used as individual detectors for matched filtering, as a single subspace detector. These detectors are subspaces made up of n eigenvectors representing the design set space, where n is some integer less than the number of templates. Detection is performed within the subspace defined by each detector. Continuous data is projected into this subspace, producing a detection statistic of between 0 and 1.
In order to divide the entire set of 637 templates used above into representative design sets, we will use a hierarchical clustering technique to create groups of similar templates based on both inter-event distance and waveform cross-correlation. These groups will then be represented by a series of subspace detectors and used in much the same way as templates are use in the matched filtering process. The results of the subspace detection will be compared to the results of the above matched filtering in order to assess the appropriateness of subspace detection in the context of a geothermal microseismicity and investigate methods of template clustering for subspace detector creation.
For this work we will continue to rely heavily on the seismic processing Python package Obspy \citep{Krischer_2015}. Matched filter detection are done using an open-source Python package called EQcorrscan \citep{Chamberlain_2014, eqcorrscan}, which allows for a scalable, multi-parallel approach to the matched filtering routine. The package’s so-called “embarassingly parallel” approach allows us to make use of the PAN computing cluster at the University of Auckland through the New Zealand eScience Infrastructure (NeSI), which cuts the processing time from months to a matter of a few days.