Anisotropy as an indicator for reservoir changes: example from the Rotokawa and Ngatamariki geothermal fields, New Zealand

We investigate the relation between geothermal ﬁeld production and fracture density and orientation in the Ngatamariki and Rotokawa geothermal ﬁelds, located in the Taupo Volcanic Zone, New Zealand using shear wave splitting (SWS). We determine the SWS parameters for 17 702 microseismic events across 38 stations spanning close to 4 yr from 2012 to 2015. We compare the strength of anisotropy to changes in ﬁeld production and injection. We also compare the orientation of the anisotropy to in situ and regional measurements of maximum horizontal stress orientation. ( S H max ). Due to the volume of unique events (approximately 160 000), shear wave phases are picked automatically. We carry out automatic SWS measurements using the Multiple Filter Automatic Splitting Technique (MFAST). The SWS measurements are interpreted in the context of stress aligned microcracks. Outside both ﬁelds and within Ngatamariki, fast polarizations align with the NE–SW regional orientation of S H max . Within Rotokawa a greater complexity is observed, with polarizations tending toward N–S. We observe increases in per cent anisotropy coinciding with the start of production/injection in Ngatamariki and then a later correlated drop in per cent anisotropy and v P / v S ratios in southern Ngatamariki as injection is shifted to the north. This relationship is consistent with pore ﬂuid pressure within the reservoir being affected by local changes in production and injection of geothermal ﬂuids causing cracks to open and close in response.


I N T RO D U C T I O N
Seismic anisotropy is a material property in which a seismic wave's velocity is dependent on its direction of propagation and the anisotropic symmetry of that system. Knowledge of the seismic anisotropy and by proxy the orientation and density of cracks in a geothermal reservoir can be useful for inclusion in numeric reservoir models and for the targeting of production and injection wells. Temporal changes in anisotropy can also indicate how the cracks and stress state within the reservoir respond to changes in field operations. These changes give information about the mechanisms underlying earthquake occurrence and the relationship between earthquakes, fluid and stress.
Recent studies have investigated anisotropy in geothermal fields associated with volcanoes (e.g. Nowacki et al. 2018) and the relation between anisotropy and fluid flow in fault zones and through hydrofracturing (e.g. Wenning et al. 2018). Yet studies of anisotropy * Now at: GFZ German Research Centre for Geosciences, Potsdam, Germany. Freie Universität Berlin, Berlin, Germany.
in operational geothermal fields are relatively rare in the refereed literature. Evans et al. (1995) measured shear wave splitting in the Geysers geothermal field in California and determined that fast polarizations were sometimes parallel to stress and other times parallel to rock fabric. They suggested that approximating permeability anisotropy with seismic anisotropy might help optimize fluid injection and steam recovery. Receiver functions were used to map a boundary between near-isotropic and highly anisotropic material that was inferred to be caused by the sudden change in compressibility of supercritical fluid in the Larderello geothermal field in Italy (Agostinetti et al. 2017). Seismic anisotropy at Coso geothermal field, California was examined with shear wave splitting (Lou & Rial 1997) to determine the fracture density and fracture orientations. Injection has also been linked directly with anisotropy at Coso with delay times being observed to increase with the onset of injection . This relationship is also observed at Krafla geothermal field in Iceland (Tang et al. 2005).
Rotokawa and Ngatamariki are two of the twenty-three geothermal systems found in the Taupo Volcanic Zone (TVZ; Chambefort et al. 2016). The two fields are located adjacent to each other Figure 1. Station locations with respect to field extent (determined from the resistivity anomaly) for Rotokawa (south/red; Risk 2000) and Ngatamariki (north/orange; Boseley et al. 2010). Black lines represent faults (Wallis et al. 2013;Litchfield et al. 2014). Grey points give earthquake locations for events that resulted in at least one high-grade shear wave splitting measurement. Seismicity is focused in northern and southern Ngatamariki as well as eastern Rotokawa as illustrated by blue density contours. Inset map (Stamen Design, under CC BY 3.0. and data are by OpenStreetMap, under ODbL) shows the location of the fields relative to the North Island of New Zealand with respect to the onshore boundaries of the whole TVZ (solid) and young TVZ (dashed) (Wilson et al. 1995).
( Fig. 1), allowing us to treat earthquakes in a combined catalogue. Both fields are actively producing electricity and together made up 27 per cent of New Zealand's geothermal power output during 2014 (Carey et al. 2015).
Power generation from geothermal fields involves the extraction of hot geothermal fluid from the reservoir (production). After energy is extracted and electricity has been generated, a portion (up to 100 per cent) of the cooled fluid is reinjected into the reservoir (injection). Injection plays an important role in reservoir management as well as providing a method for wastewater disposal (Stefansson 1997). Production and injection influence anisotropy by altering the pore-fluid pressure and effective stress field, causing cracks to open/close and potentially form. The microseismic events that we use as a seismic source are also largely induced by injection (Eberhart-Phillips & Oppenheimer 1984;Sherburn et al. 2015).
With ongoing seismic monitoring across both fields, we have an earthquake catalogue spanning close to 4 yr from April 2012 to November 2015. 17,702 earthquakes, with their S-wave arrivals picked, give 49,967 station-event pairs across 38 stations that we analyse for shear wave splitting. To determine how the shear wave splitting fast polarizations vary across the fields, we apply a spatial averaging technique . We also inspect the time variation in shear wave splitting delay times with relation to changes in production and injection within the fields. Wellhead locations for injection (red/diamond) and production (blue/square) wells in Ngatamariki. The field extent (determined from the resistivity anomaly) is outlined in orange (Boseley et al. 2010).

Field operations and seismicity
The majority of the seismicity in both fields is induced by stress changes from the injection of colder geothermal fluid through a combination of rock cooling/contraction and pressure changes (Eberhart-Phillips & Oppenheimer 1984;Sherburn et al. 2015;Hopp et al. 2019).
Injection and production in Ngatamariki largely occurs at depths between 2-3 km (with respect to sea level) and is centred between the majority of the seismicity, which ranges in depth from 1.5 to 3 km (Fig. 1). Production wells are located in the central-southern portion of the field (Fig. 2). Injection takes place in two separate locations, in the north and in the south of the field. Seismicity is divided into a northern cluster and a southern cluster associated with the respective injection areas.
Within Rotokawa, production takes places to the northwest of the field (∼2.5 km depth) and, divided by the Central Field Fault, deeper (up to 3 km depth) injection takes place in the southeast Fig. 3). The single major cluster of seismicity is associated with the southeast injection wells and is focused from 1.5-2.5 km depth (Fig. 1).

Shear wave splitting and aligned cracks
Shear wave splitting occurs in an anisotropic medium when the component polarized in one particular orientation travels faster than its orthogonal component. Shear wave splitting is characterized by two variables: the polarization of the leading wave (φ) and the delay time (δt) between the fast and orthogonal polarizations (Crampin 1978;Savage 1999) Shear wave splitting in the upper crust is often caused by aligned microcracks parallel to the orientation of maximum horizontal stress Downloaded from https://academic.oup.com/gji/article-abstract/220/1/1/5561446 by Lawrence Berkeley Lab user on 03 March 2020 . Wellhead locations for injection (red/diamond) and production (blue/square) wells in Rotokawa. The field extent (determined from the resistivity anomaly) is outlined in red (Risk 2000). (S Hmax ;Crampin & Peacock 2008). The underlying principle is that cracks not aligned with maximum horizontal stress are closed and any fluid they contained is forced into the cracks aligned with S Hmax . We include pore-space in our definition of microcracks.
Cracks and pores also affect the ratio of P-and S-wave velocities, or the v P /v S ratios, so comparing splitting with v P /v S ratios may help to distinguish different processes (e.g. Lucente et al. 2010). But the relation between anisotropy and v P /v S ratio is complicated by the different effects of aligned cracks, pore volume, and also the type of the material in the cracks. Shear wave splitting depends only on S waves, which respond similarly to liquid or gas-filled cracks. Large amounts of liquid-filled pores or cracks increases v P /v S ratios (generally values >2 suggest high amounts of liquid). However, P-wave velocities are lowered dramatically with gas-filled cracks because of the high compressibility of gas, leading to low v P /v S ratios (generally on the order of 1.6 or lower; e.g. Chatterjee et al. 1985). The sign of the correlation between v P /v S and per cent anisotropy may switch depending on whether the cracks are saturated or unsaturated (Unglert et al. 2011). Saturated cracks result in a positive correlation while unsaturated cracks result in a negative correlation. For saturated cracks, v S reduces faster with increasing crack density than v P resulting in an increase in v P /v S . For unsaturated cracks, v P decreases faster with increasing crack density than v S resulting in a decrease in v P /v S . When cracks become saturated v P increases (with smaller increase in v S ) which results in an increase in v P /v S . Alternatively, a negative relation between delay time and v P /v S ratio could be caused by changing shapes of the cracks-higher pressures can narrow the cracks, leaving less pore space and so lowering v P /v S ratios, but increasing anisotropy because of more highly aligned cracks (Lucente et al. 2010;Zal et al. 2018).

Previous results in the Taupo Volcanic Zone
S Hmax orientations have been determined across the Taupo Volcanic Zone (TVZ) from focal mechanisms by Townend et al. (2012). They report that the average orientation within the TVZ is NE-SW or ENE-WSW in a prominently normal faulting regime (Hurst et al. 2002). These orientations provide an expected average shear wave splitting fast polarization without any prior knowledge from the fields themselves.
Three wells within Rotokawa have also had their maximum horizontal stress orientations measured in 2010 (McNamara et al. 2015) and orientations around NE-SW are reported. These provide an important comparison for the shear wave splitting polarizations that we measure. Fracture width and spacing distribution within these wells has also been studied (Massiot et al. 2017).
Shear wave splitting in the TVZ has been determined for crustal earthquakes in several studies. In a broad study of the North Island (Audoine et al. 2004), shallow events measured at high frequencies generally had fast polarizations oriented with the strike of the TVZ, as expected for extensional cracks opening parallel to the TVZ and consistent with the results of the stress inversions. Near Ruapehu volcano fast polarizations were mixed, possibly due to structurally varying anisotropy or to varying stresses near the tip of the TVZ (e.g. Johnson 2011). A transect across the TVZ north of Lake Taupo, near this study region, used deeper earthquakes to find rift-perpendicular fast polarizations within the TVZ (Morley et al. 2006). P-wave anisotropy determined from tomography yields fast orientations in this region at shallow depths (4 km) that have rapid spatial variations (Eberhart-Phillips & Reyners 2009).

Station orientation
Most of the stations were deployed on the surface using a magnetic compass to determine their orientations. However, three borehole stations within Ngatamariki (NS12, NS13 and NS14) have not previously had their orientations calibrated. Thus any shear wave splitting polarization measured on these stations will have a fixed offset added that is equal to the misorientation. To account for this we determine the correct orientation of the stations and then apply the correction to the fast polarizations post processing.
We orient the stations using the Rayleigh-wave method of Stachnik et al. (2012), which was applied in orienting ocean bottom seismographs in Cascadia (e.g. Lodewyk & Sumy 2015). The method relies on the particle motion of Rayleigh waves, in which the vertical and radial components are out of the phase by π /2, so that the Hilbert transform of the vertical component should be nearly identical to the radial component. Trial values of the orientation of the seismometers are used to calculate the radial components and to measure their similarity with the vertical components. The best orientation for any particular earthquake-station pair is the one that maximizes the correlation coefficient between the radial and vertical components. We take all earthquakes in the 2014/2015 USGS catalogue, with magnitude greater than 6 and a depth of less than 100 km, between 30 • and 180 • epicentral distance from Rotokawa and Ngatamariki. Waveforms are 10 per cent cosine tapered with a bandpass filter from 0.02 to 0.04 Hz applied. Results are summarized in Table 1 including standard error calculated with a bootstrap method similar to that described in the Supporting Information (Section S1) for mean fast polarization but with the angles not doubled.

Earthquake catalogue
Our catalogue is made up of earthquakes detected by GNS Science (e.g. Sherburn & Bourguignon 2015) and supplemented with earthquakes detected using the matched filter method (Hopp et al. 2016;Chamberlain et al. 2017;Hopp et al. 2019). The catalogue spans almost 4 yr from March 2012 to September 2015. Earthquakes in our catalogue have P-wave picks made either automatically, for the GNS Science events, in the seis-comp3 software package or, for matched filter events, based on the pick of a detection's template (Hopp et al. 2019; Section 2.3). We automatically pick S-wave arrivals (Section 2.4). The ratio of the number of microseismic events located in Rotokawa compared to Ngatamariki is approximately 2:1. Station completeness is inspected in the Supporting Information (Section S2). At Ngatamariki, several stations were commissioned at the same time as production began in early 2013, and station coverage also changed in early 2015. In Rotokawa station coverage change in late 2012 and early 2013, but was fairly consistent after that. All earthquakes are relocated with the program NonLinLoc (Lomax et al. 2009) using the VELEST (Kissling et al. 1994) minimum 1-D velocity model determined for the fields (Sewell et al. 2017; Table 2). For depths greater than 3 km (deeper than the majority of the seismicity), velocities are averaged from the GNS 3D model for Rotokawa (Rawlinson 2011,p. 77).

Matched filter detection
The matched filter (also referred to as template-matching) method of detecting earthquakes relies on waveform cross-correlation of a known (template) event with continuous seismic records (Shelly et al. 2007). When a detection statistic (the sum of cross-correlation coefficients for each channel) reaches a threshold value (here we use eight times the mean absolute deviation of the statistic) a detection is made. The detections within our catalogue are made using the EQcorrscan Python package (Chamberlain et al. 2014).
The matched filter method is particularly suited to detecting repeating signals in continuous data with high levels of noise (Gibbons & Ringdal 2006). Therefore, it is ideal for use in geothermal areas with high levels of natural (e.g. Clacy 1968) and anthropogenic noise related to fluid injection and production. This method proved to be effective at Rotokawa and Ngatamariki with 6384 template events resulting in 154 635 additional detections.
Most events detected using the matched filter method are of poor quality (hence their failure to be detected by more conventional methods). While approximately 1 000 000 station-event pairs were detected this way, only 2.1 per cent of these resulted in a successful S-wave pick. Due to the large volume, adding the detected events approximately doubled the catalogue of station-event pairs that could be analysed for shear wave splitting. Approximately an equal number of template and detected events resulted in a high-grade shear wave splitting measurement. Stacking detected events from a single family (detected from one template) could be a future path of research particularly to measure the time variation within a single family.

Automatic shear wave picking (spickerC)
Shear wave splitting measurements require the arrival time of a shear wave phase to be determined. The volume of events in the catalogue make manual picking time prohibitive so instead an automatic phase picker is used, spickerC of Castellazzi et al. (2015).
The automatic shear wave picking program (S-picker) was originally written by Diehl et al. (2009). The version we use here was modified from version 1.4.0 (for use with local earthquakes in New Zealand) by Castellazzi et al. (2015). This version of S-picker (1.4.0.a) is known as spickerC. spickerC was originally used by Castellazzi et al. (2015) in conjunction with MFAST (Section 2.5) to make automatic shear wave picking and splitting measurements at Ruapehu Volcano in New Zealand.
spickerC uses a combination of three different detection and picking methods: short-term average versus long-term average ratio (STA/LTA) as described by Allen (1978), polarization detection based on the approach of Cichowicz (1993), and an autoregressive picker using the Akaike Information Criterion (AR-AIC) as described by Leonard & Kennett (1999). The average pick is a weighted average of the three methods with a final grade based on the error estimate and the signal-to-noise ratio around the pick (Castellazzi et al. 2015). We keep the picks that are assigned with the highest grade (class0). The grading scheme and processing parameters for spickerC are given in the Supporting Information (Section S3). To account for shear wave splitting we use the earliest AR-AIC pick instead of the average pick. Supporting Information Fig. S3 shows an example of the different S-wave picks for a single event.

Automatic shear wave splitting (MFAST)
We use the Multiple Filter Automatic Splitting Technique (MFAST; Teanby et al. 2004;Savage et al. 2010;Wessel 2010). MFAST conducts shear wave splitting measurements following the method of Silver & Chan (1991), with corrections for error bars from Walsh et al. (2013), along with the clustering analysis of splitting measurements of Teanby et al. (2004).
For this paper we use the original MFAST software package rewritten into the R software environment (R Core Team 2013). Based on MFAST version 2.2, MFASTR takes advantage of R's capabilities, most notably in parallel computing, to greatly increase Downloaded from https://academic.oup.com/gji/article-abstract/220/1/1/5561446 by Lawrence Berkeley Lab user on 03 March 2020 processing speed and to provide added flexibility. While MFASTR mirrors MFAST in most respects, it does have some key processing differences which are outlined in the Supporting Information (Section S4) along with the main processing steps of MFAST.
To grade shear wave splitting measurements we follow a scheme similar to Castellazzi et al. (2015) where dependable measurements should be consistent over a number of different filters. Measurements that are graded F1, F2 and F3 are termed 'high grade' and are given weighting of 1, 2 and 3, respectively. The details of grading and weighting schemes are given in the Supporting Information (Section S4).

Events with shallow propagation angles
Ideally, the shear wave splitting measurements would only contain events with propagation angles (with respect to the anisotropic symmetry axis) such that their delay time is a maximum and their fast polarizations are subparallel to the strike of the cracks (i.e. Band 2 in Peacock et al. 1988). MFASTR, in addition to angles of incidence, records the ray parameter for each event. Thus we are able to use the ray parameter and the velocity model to remove events whose propagation angles make the measurement unreliable. This also works to avoid phase interference at the surface by removing events outside of the shear wave window.
For determining ray parameters and incidence angles, we use the same velocity model that the events are located with (Section 2.2).
For our catalogue of shear wave splitting measurements we remove all events that did not spend at least 17 per cent of their path with a propagation angle less than 35 • . This is estimated by firstly finding the depth at which propagation angle goes below 35 • . The straight line distance from this point to the hypocentre is then calculated and given as a percentage of the total straight line path. If it exceeds 17 per cent then the event is excluded.
We choose a 17 per cent cut-off by making a two-layer approximation with a maximum velocity of 2 km s −1 and per cent anisotropy of 5 per cent in the final layer and 2.5 per cent in the lower layer. Numerical experiments show that if 17 per cent or more of the total path is in the final layer, then we measure a fast polarization that corresponds to the crack orientation in the top layer.
By taking the measurements that fall within the 17 per cent cutoff we are able to make an estimate of strike of the crack plane that the measurements are sampling from (by calculating their mean fast polarization). The anisotropic symmetry of aligned vertical cracks means that measurements with backazimuths close to the strike of the cracks have delay times and polarizations that are insensitive to their propagation angle. To account for this we also include measurements that have backazimuths close to the mean fast polarization of the station they are measured on (calculated from events with sufficiently steep propagation angles). To achieve this we take the difference between the backazimuth and the fast polarization (φ diff ) along with the take-off angle for that particular event (θ take-off ). By calculating and taking an arbitrary cut-off angle of 15 • (θ cut ). An event meets this criteria if In practice this means that any event with take-off angles less than 15 • meet the criteria regardless of their backazimuth. Events with greater take-off angles must have backazimuths closer to the mean crack orientation. The limit for very shallow take-off angles is that backazimuths must be within 15 • of the mean crack orientation.

Multimodal fast polarization distributions (Kuiper's Test)
We use the Rayleigh test for circular uniformity to identify single stations that do not have a dominant circular mean (i.e. a uniform distribution of fast polarizations). The Rayleigh test is further described in the Supporting Information (Section S1).
In the case that a set of fast polarizations (usually measured on a single station) fail the Rayleigh test, Kuiper's one sample test is applied. Unlike the Rayleigh test, Kuiper's test examines multimodal alternatives to uniformity. In essence, a sample that fails the Rayleigh test but passes Kuiper's test does not have a dominant mean polarization, but may have several modes. If Kuiper's test is passed then the measurements are used in spatial averaging. If Kuiper's test fails we treat this as evidence that the polarizations carry little or no information. It is conceivable that the fast polarizations are changing with backazimuth in such a manner that together they appear uniform. To rule out this possibility we apply a weighted correlation test between the (doubled) fast polarizations and the backazimuths. If there is little or no evidence (p-value ≥ 0.05) that the correlation is finite then we discard the corresponding measurements as they do not carry any information. To apply this correlation we use the method of Agostinelli & Lund (2013) which is the circular equivalent of Pearson's product-moment correlation.

Spatial averaging (TESSA)
Spatial averaging of shear wave splitting fast polarizations gives a quick and generally good way of digesting thousands of measurements at once to understand the average trends of the data set. For this approach we use the TESSA (Tomography Estimation and Shear-wave-splitting Spatial Average) package (Johnson 2011). TESSA includes codes for computing a 2-D tomography estimate for delay time, however we instead opt for an approach that focuses more on temporal variation rather than spatial variation of delay times (Section 2.8). In this section we summarize the method of Johnson (2011);Johnson et al. (2011) for the spatial averaging of fast polarizations.
We apply an artificial weighting to TESSA by repeating measurements (i.e. F2 and F3 grade measurements are repeated twice and three times respectively while F1 grade measurements are only included once).
Due to our much smaller study area our parameters differ from those used by Johnson et al. (2011). To reflect our study area being roughly a fifth of the size we adjust most parameters by a factor of 5. Original parameters can be found in Johnson et al. (2011) and our parameters are summarized in Table 3.
For the full catalogue, we grid our study area using the quad-tree method (Townend & Zoback 2001). The gridding is performed with respect to the maximum and minimum allowable number of rays per cell as well as the minimum cell size (Table 3). The quad-tree algorithm works in an iterative process. At each iteration the number of ray paths in a cell is counted and if it is above the maximum the cell is divided further and if it is below the minimum the cell is excluded (Johnson 2011, Chapter 4). This process continues until either the number of rays falls between the minimum and maximum or the cell size reaches the minimum allowed (Johnson 2011, Chapter 4).
When making temporal comparisons with our data set we use a regular grid so that the grid is the same for all the time intervals. The size of the regular grid is based solely on the settings for minimum cell size 'minsize' ( Table 3). Rays that pass through each grid cell are weighted and then averaged using the correct considerations for axial quantities (Supporting Information Section S4.3). We weight measurements by distance (d) along ray path by 1/d 2 as the fast polarization will be affected more by anisotropy later in its path (Rümpker & Silver 2000;Nistala & McMechan 2005). The 1/d 2 weighting scheme gave the lowest scatter of those tested by Johnson et al. (2011).

Determining temporal variation in per cent anisotropy and v P /v S
To determine how delay times and per cent anisotropy vary through time we divide the study period into week long windows, determine the per cent anisotropy for each event in the window, and then calculate a mean value of per cent anisotropy for that week. Events are weighted by their grade (Supporting Information Section S4.3) and events with shallow propagation angle are removed (Section 2.6). The one week window size is chosen so a sufficient number of measurements are averaged (greater than 20) in each window and so that transient earthquake swarms do not dominate.
To aid in the interpretation of the raw time-series we also perform a cumulative sum and then remove the trend line (by fitting a linear model). In effect we are observing the deviations in weekly per cent anisotropy from the mean. A negative slope corresponds to a period where values are less than the mean, a positive slope where values are greater than the mean, and a flat slope where the values are at the mean (and the flatter the slope the closer the values are to the mean). The cumulative sum also acts like an integrating (low-pass) filter and so attenuates short-period variation. This cumulative sum does not have any physical meaning for per cent anisotropy and serves solely as a method of filtering the data. We refer to this measure of cumulative per cent anisotropy as k cum .
These steps are also followed with v P /v S ratios (P-wave velocity to S-wave velocity ratios). v P /v S can be calculated for all events with an S-wave pick by dividing the S-wave traveltime by the P-wave traveltime. Origin times from the event locations, along with the Pwave (seiscomp3 and eqcorrscan) and S-wave (spickerC) picks, are used to determine the traveltimes. We determine v P /v S for all events with a class0 S-wave pick grade (Section 2.4). v P /v S is sensitive to uncertainties in origin time, P-wave pick time, and S-wave pick time. In the Supporting Information, we manually pick a subset of events to determine if the automatic picking is stable through time (Supporting Information Section S5). Our method of averaging over many events and, particularly, performing the cumulative sum leaves the resulting time-series relatively insensitive to large uncertainties in the origin and pick times.

Single station analysis
Before presenting results from spatial averaging, it is important to understand trends for the stations that the measurements originate from. Fast polarizations at stations surrounding Rotokawa and Ngatamariki are plotted in Figs 4 and 5 respectively. We include blue bars to indicate the mean for stations that pass the Rayleigh test (pvalue < 0.05) as there is evidence for a unimodal departure from uniformity and thus a single mean axis carries information. Stations that fail the Rayleigh test but pass Kuiper's test (Section 2.6.1) may carry visual information so are shown just by their rose diagrams (e.g. RT19). Stations that fail both tests carry no visual information (their rose diagrams would be a circle) and thus are represented simply by their location.
Supporting Information Table S3 summarizes the (weighted) shear wave splitting parameters for stations whose fast polarizations pass the Rayleigh test. The full results for each station-event pair are in the Github repository described in the Acknowledgements. Supporting Information Table S4 summarizes stations that fail the Rayleigh test. No fast polarization is quoted for these stations as it carries no information. Events with shallow propagation angle are removed based on their ray parameter (Section 2.6). During the study period two stations were disbanded and moved to the nearest fence-line after being disturbed. For measurements reported by station, we merge NS08 into NS18 and NS09 into NS16. The number of measurements we have before the stations are moved is low enough that they have little impact. The true station locations are persevered in spatial averaging (Section 3.2).
Of the 1.2 million station-event pairs in our catalogue, we successfully pick the S-wave arrival for 50 000 (4 per cent) of these. Of those with an S-wave pick, we get high grade shear wave splitting measurements from 13 000 (26 per cent). So approximately 1 per cent of the events located on a particular station result in a high grade shear wave splitting measurement. Castellazzi et al. (2015) recorded a much higher rate of S-wave picks (28 per cent). This dramatic difference is due to the matched filter detected events in our catalogue. The template events have a comparable picking rate of 34 per cent. The largest difference between the two studies is in the percentage of S-wave arrivals resulting in a high-grade measurement. For the Ruapehu study Castellazzi et al. (2015) find 60 per cent of S-wave arrivals led to high-grade measurements versus 26 per cent herein. This is probably due to differences in study region or higher quality S arrivals in the Ruapehu study. We focus on a small area with many low magnitude and thus lower quality earthquakes. In contrast, Castellazzi et al. (2015) use larger earthquakes (magnitude >1.5) in a much larger study area.

Borehole stations
The three borehole stations (Section 2.1) report mean fast polarizations close to what we expect based on nearby stations (Supporting Information Table S3, Fig. 5). Particularly, NS13 (43 • ) and NS14 (52 • ) have mean fast polarizations close to the prevailing stress orientation of NE-SW (Townend et al. 2012) as well as having means close to those of nearby stations. This gives us confidence that the station orientations that we have determined are accurate. NS12's 5 • mean fast polarization is further from the expected 45 • . However, its mean fast polarization is still similar to other surface stations in the area (e.g. NS18, NS02, RT22 and NS07) so we consider its orientation to be correct. . Fast polarization rose diagrams for Rotokawa stations for events that satisfy incidence angle criteria (Section 2.6). Mean weighted fast polarizations for stations that pass the Rayleigh test (p-value < 0.05) are marked with blue lines. Stations with uniform distributions are indicated with blue triangles and those that have multimodal departures from uniformity (pass Kuiper's test) just include their rose diagram. The field extent (determined from the resistivity anomaly) for Rotokawa is outlined in red (Risk 2000).

Figure 5.
Fast polarization rose diagrams for Ngatamariki stations for events that satisfy incidence angle criteria (Section 2.6). Mean weighted fast polarizations for stations that pass the Rayleigh test (p-value < 0.05) are marked with blue lines. Stations with uniform distributions are indicated with blue triangles and those that have multimodal departures from uniformity (pass Kuiper's test) just include their rose diagram. Borehole stations are marked with * . The field extent (determined from the resistivity anomaly) for Ngatamariki is outlined in orange (Boseley et al. 2010).

Per cent anisotropy
We can also determine the mean per cent anisotropy (Savage 1999) for each station as well as an average for Rotokawa and Ngatamariki stations separately. β max , β min andβ correspond to the maximum, minimum and mean S-wave velocities, respectively: β is estimated from the straight line path distance (d), traveltime (t) and the delay time (δt): Supporting Information Table S5 shows the per cent anisotropy and differential shear wave anisotropy (Crampin 1999) for all stations whose events pass the Rayleigh test (Supporting Information Section S1). Mean values are determined by assuming a homogeneous medium (i.e. straight line paths with velocities determined from traveltime and delay time).
For means across Rotokawa (RT) stations we find that k = 4.00%.
For Ngatamariki stations (NS) the mean is approximately 0.5 per cent less than that of Rotokawa stations: These values represent the average in the top 3 km or so of the crust and correspond well with other crustal studies (Crampin & Peacock 2008). It is likely an underestimate of the near surface as many shallow events are removed due to their shallow propagation angles. Values for each station (Supporting Information Table S5) fall into the range expected for an intact rock (Crampin 1999).
The mean delay time of the full data set is 0.12 s with a standard deviation of 0.073 s.
Average values of v P /v S ratios are around 1.6, consistent with the large amount of gas-filled cracks and pores in geothermal areas (Chatterjee et al. 1985;Unglert et al. 2011).

Spatial averaging
To display the spatial averaging results from TESSA we employ a slightly different approach to that of Johnson et al. (2011). Instead of removing the rose diagrams from grid cells that have a standard error and standard deviation greater than some value, we apply the weighted Rayleigh test (Supporting Information Section S1) using the weights that TESSA outputs. If that grid cell fails the Rayleigh test then it is not displayed.

Full catalogue
We get the best resolution by applying the spatial averaging to the entire catalogue (illustrated in Fig. 6a). By using the quadtree gridding scheme (Section 2.7), we can compare regions with dense seismicity to surrounding regions with less events. We are effectively averaging through time and the results represents the most dominant S Hmax orientations throughout the study period.
Within Ngatamariki, northern Rotokawa, and to the west of both fields we observe mean fast polarizations close to the expected NE-SW S Hmax orientation (Townend et al. 2012).
There is also a rotation of polarizations in northern Rotokawa into the Aratiatia fault zone (Fig. 7). However, this is only apparent on one station (RT13; Fig. 4) and so should be treated with caution as it is possible that it is the product of RT13 being misorientated or a local effect.
Within Rotokawa (as well as to the east) and to some extent central Ngatamariki there is an N-S trend to the mean fast polarizations. McNamara et al. (2015) observe a subset of cracks striking N-S within their studied wells (Section 3.3) so it is possible that there are areas where these N-S cracks are dominant. McNamara et al. (2015) do not note any scale difference between the N-S cracks and other crack families. While these N-S orientations are also somewhat unexpected, they are present on more than one station (e.g. RT12, RT22, RT09 in Fig. 4), so we are more confident that they are real.

Temporal variation in fast polarization
For investigation of temporal variation in fast polarization using spatial averaging, we use a regular grid rather than a quadtree grid (Section 2.7). By using the regular gridding scheme, the grid remains the same between time periods and thus makes the comparison easier and fairer .
We break our time period into three portions that divide the study period into (b) before production began in Ngatamariki (March 2013), (c) the time period between the beginning of production and the commissioning of NM12 production well (October 2014), and (d) the period from the commission of NM12 to the end of the catalogue. These are plotted in panels (b)-(d) of Fig. 6.
Overall variation between years appears minimal if we ignore the large rotation into the Aratiatia fault zone within northern Rotokawa from station RT13. N-S polarizations remain ubiquitous in eastern Rotokawa.

Relation to faults
The TVZ is largely a normal faulting regime (Hurst et al. 2002) thus we expect the faults to generally strike in the same orientation as that of S Hmax . This means we do not expect to be able to distinguish well between anisotropy resulting from fault parallel shear (Savage et al. 1990;Zhang & Schwartz 1994;Evans et al. 1995) and anisotropy originating from fluid filled microcracks (from just orientation differences).
In northern Rotokawa (Fig. 7), within the Aratiatia fault zone, (for the full catalogue) we see fast polarizations that lie almost parallel to the strike of the faults (ignoring measurements from RT13 that are heavily rotated into the E-W (Fig. 4)). Within central Rotokawa there is a tendency for mean fast polarizations to be close to the strike of the major faults. Mean fast polarizations rotate towards the north and closer to the strike of the production field fault (PFF, mean: 20 • ), central field fault (CFF, mean: 30 • ), and injection field fault (IFF, mean: 30 • ) (Wallis et al. 2013; Fig. 7).
Temporally, overall these relationships with the faults are maintained (Fig. 6). Particularly, the production field Fault (PFF) has mean fast polarizations close to its strike after 2012. Again, ignoring RT13, fast polarizations surrounding the Aratiatia fault zone (northern Rotokawa) are also close the strike of the nearby faults.

Relation to in situ stress measurements
Three wells within Rotokawa have had their horizontal stress orientations determined by McNamara et al. (2015; Fig. 7). These measurements are of great importance to this study as they provide a direct in situ measurement of S Hmax that we can compare against. The logs were made in 2010, approximately 2 yr before our earthquake catalogue begins. Fig. 7 illustrates the mean S Hmax orientations for the three wells determined from drilling induced tensile fractures, petal centreline fractures, and borehole breakout.
Examination of Fig. 7 shows that the spatially averaged fast polarization in the cells adjacent (within 300 m) to RK32 are close to the orientation measured by McNamara et al. (2015). RK18L2 differs more strongly, however this is interpreted as due to activity on the Production Field Fault that it intersects (McNamara et al. 2015). All three wells have means that are statistically different from the mean fast polarization of the nearest spatial averaging cell (p-values < 0.05). In terms of absolute values all are within 20 • (19 • , 13 • and 7 • for RK30L1, RK32 and RK18L2, respectively). Given that these horizontal stress orientations predate our catalogue and the absolute differences are small, the in situ measurements are compatible with what we observe in the shear wave splitting fast polarizations. Notably, between RK32 and RK18L2 there is a 20 • rotation of the mean horizontal stress orientation while the wellhead locations are only separated by approximately half a kilometre (a significant rotation over a comparatively small distance; McNamara et al. 2015). Most of the variation we observe in the Rotokawa fast polarizations is on this scale and so may be reflecting real changes in the underlying stress orientations rather than uncertainty in the shear wave splitting measurements.

Temporal variation in per cent anisotropy for Ngatamariki
Production began at Ngatamariki in early March of 2013. This gives us a rare opportunity to assess how the anisotropy within the Ngatamariki reservoir reacts to the start of the field being stimulated. Seismicity within Ngatamariki is divided into a northern and Downloaded from https://academic.oup.com/gji/article-abstract/220/1/1/5561446 by Lawrence Berkeley Lab user on 03 March 2020 Figure 6. Spatially averaged fast polarizations for all high-grade measurements weighted by grade. Grid cells that pass the Rayleigh test (p-value < 0.05) are plotted. Black lines represent faults (Wallis et al. 2013;Litchfield et al. 2014). The field extents (determined from the resistivity anomaly) for Rotokawa to the south and Ngatamariki to the north are plotted in red and yellow, respectively. Map panels display results for (a) the full catalogue, (b) before production began in Ngatamariki (March 2013), (c) the time period between the beginning of production and the commissioning of the NM12 production well (October 2014) and (d) the period from the commission of NM12 to the end of the catalogue. a southern cluster (Fig. 1). Injection occurs in the north and the south while production occurs in the southern to central portion of the field. Given this obvious north-south division we divide our catalogue into two portions such that ray paths cross either the northern or the southern portion of the field. Shear wave splitting has been observed to be sensitive to fluid pressure in a geothermal setting (e.g. Tang et al. 2005) and so we expect some relationship between delay times and injection/production.
To rule out the migration of seismicity as a cause for changes in anisotropy, in the Supporting Information, we also limit the event locations to the southern cluster and the northern cluster for southern and northern Ngatamariki respectively (Supporting Information Section S6) and find that the relationships are more pronounced or remain the same.

Southern Ngatamariki
All of the production and a significant proportion of the injection occur in the southern-central portion of Ngatamariki (Fig. 2). Most of the seismicity is also confined to the south (Fig. 1). For analysis of the southern portion of the field, we take all events on the stations located on or near the southern locus of seismicity (NS11, NS06 and NS14; Fig. 5). For stations to the south of the field (RT22, NS07, NS10, NS15 and NS16) we take all events north of −38.58 • latitude. And finally, for stations near the centre of the field (NS02, NS05, NS12 and NS13) we take all events south of −38.55 • latitude. We include events from Rotokawa because shear wave splitting is generally more sensitive to anisotropy later in the station-event path (Rümpker & Silver 2000;Nistala & McMechan 2005) and in this section we are investigating relative changes rather than absolute. Fig. 8 shows the temporal results for southern Ngatamariki. After production begins we see a higher than average anisotropy followed by a stable period near the mean during most of 2014. In late 2014 an increase in injection and a larger magnitude increase in production precedes a correlated drop in anisotropy and v P /v S ratios.

Northern Ngatamariki
For northern Ngatamariki we take all events on stations north of the northern cluster of seismicity (NS01, NS03, NS18 and WPRZ; Fig. 5). On the stations NS02, NS12, NS13, NS04 and NS05 we take all events north of −38.55 • latitude. Fig. 9 shows the temporal results for northern Ngatamariki. As with southern Ngatamariki, initially there is higher than average anisotropy in response to the beginning of injection and production. No production occurs in the northern area and, beginning in early 2015, uncertainty in anisotropy increases. Unlike in the south, v P /v S ratios remain relatively constant throughout the study period.

Temporal variation in per cent anisotropy for Rotokawa
For the analysis of the temporal variation in per cent anisotropy within Rotokawa we take all events occurring south of Ngatamariki (south of −38.58 • latitude) measured on all Rotokawa stations (RT) and NS10 (the southernmost Ngatamariki station). Station locations are shown in Fig. 4 and wellhead locations in Fig. 3. Fig. 10 shows the temporal results for Rotokawa. Stable injection volume is associated with anisotropy that remains relatively constant throughout the study period. v P /v S shows larger changes that may be anticorrelated with the anisotropy, suggesting differing mechanisms from Ngatamariki.

Ngatamariki temporal variation
A cumulative plot of a time-series act as a low-pass (integrating) filter. Shear wave splitting delay times (and hence per cent anisotropy) are highly scattered (Crampin & Peacock 2008). Filtering out the dominating high-frequency variation allows longer term trends to be identified.

Southern Ngatamariki
What is immediately apparent from the plot of the temporal variation in anisotropy and production (Fig. 8) is the rapid increase in k cum (cumulative per cent anisotropy) several months after production begins. This increase levels offs and drops slightly as production and injection in the south of Ngatamariki stabilizes. Anisotropy increased when injection began at both Coso  and at Krafla (Tang et al. 2005), suggesting that it is a general feature at the start of geothermal production. The beginning of injection and its related induced seismicity may be opening new or pre-existing cracks that are favourably oriented to the stress field, therefore increasing the anisotropy by increasing the number of open cracks.
In late 2014 there is a large increase in production and a smaller increase in injection after the commissioning of a new production well (NM12). Overall production/injection is being increased, however the northern injectors take up some of the increased injection. This precipitates a rapid drop in k cum which corresponds to lower than average per cent anisotropy. This drop in anisotropy, beginning late 2014 to the end of the catalogue, may be due to the additional fluid being removed from the southern portion of the reservoir. This fluid is then being injected into the north where the connection to Date Figure 8. Southern Ngatamariki: weekly mean per cent anisotropy (k, black/solid) and weekly mean v P /v S ratios (v P /v S , red/dashed) are plotted in both their cumulative (acting as a low-pass filter) and unfiltered forms. Daily injection/production (tonnes) is plotted in blue/dot-dashed and green/dashed, respectively (with a 7-d median filter applied). 95 per cent confidence intervals are indicated by dotted lines and the vertical blue/dashed line marks the beginning of production. Labels on the x-axis mark the start of that year. An increase in production and injection in late 2014 precedes a decrease in anisotropy.
the south is not strong (Buscarlet et al. 2015), causing pressures within the southern part of the reservoir to drop and consequently close cracks. Given that we assume that the anisotropy we measure (or at least the component of it that is changing) comes from aligned fluid filled cracks, v P /v S ratios can provide another measure of the crack properties in the southern portion of the reservoir (Fig. 8). Interestingly, cumulative v P /v S ratios do not respond to the beginning of production unlike the anisotropy (the rapid increase in cumulative v P /v S at the beginning of the catalogue is a result of sparse data causing v P /v S to be larger than the mean of the full time period). Invariant v P /v S suggests that the cracks which are responding to the beginning of production are only the stress aligned cracks that we measure with shear wave splitting (while remaining invisible to v P /v S which is calculated from the traveltime of the fast split wave only).
However, in late 2014 (when there is a spike in production) it appears as if v P /v S becomes correlated with anisotropy. They follow a similar downwards trend (marked with initially sharp increase in v P /v S ) which corresponds to both v P /v S and per cent anisotropy being lower than their respective means during this period. This behaviour is consistent with increased production causing the proportion of cracks filled with steam to increase through boiling (decreasing v P /v S ). Pressure drops associated with increased mass extraction would also close cracks, causing a correlated drop between v P /v S and per cent anisotropy. A new equilibrium between pressure and crack density would be reached with a higher fraction of steam filled cracks in the reservoir and v P /v S and per cent anisotropy would stabilize. This correlation indicates that both the stress aligned cracks and the randomly aligned cracks are both responding to the additional removal of fluid from the southern portion of the reservoir.
This leads to the hypothesis that when production/injection begins in a previously unexploited geothermal reservoir the cracks which are already aligned in the stress regime are opened at the expense of unfavourably oriented cracks. Given that this spike in per cent anisotropy as production begins is relatively short lived it appears as if the reservoir reaches some equilibrium as production and injection stabilize and the cracks that were forced open by the sudden beginning of stimulation return close to pre-stimulation levels.

Northern Ngatamariki
Within northern Ngatamariki only injection is occurring and thus will be the main operational factor that will affect the cracks present.
Initially we see a similar trend to southern Ngatamariki where the beginning of production appears to cause an increase in k cum , however the change is not as strong as in the south (Fig. 9). However, in early 2015 there is a further increase in the k cum . The error bars during this latter period are large, but if the trend is real, then because this section of the reservoir is not directly being produced from, our hypothesis here is that the mechanism by which the k cum continues to increase is the same as we see at the beginning of production in both the northern and southern portions of the reservoir. That is, the stress aligned cracks being held open by the injected volume.
Because there is no production in the northern area of the reservoir the cracks remain open. This idea is corroborated by the v P /v S ratios (Fig. 9). When compared to southern Ngatamariki (Fig. 8), the v P /v S remain remarkably constant (not deviating far from its mean) even through the increase in injection and subsequent increase in anisotropy. v P /v S and anisotropy remain largely uncorrelated, which suggests that the cracks that are responding to changes in injection in northern Ngatamariki are the stress aligned cracks that we measure with shear wave splitting (Section 4.3). Since only injection is occurring in northern Ngatamariki, the anisotropy is being driven by the volume of water being injected. In contrast, in the southern portion of the field the anisotropy is initially being driven by the injection, however (at what appears to be a certain threshold) the net removal of mass results in a drop in pore-fluid pressure followed by the cracks closing and thus a decrease in anisotropy that correlates with the v P /v S ratios.

Rotokawa temporal variation
Within Rotokawa we do not have any major changes in production and injection for the duration of our catalogue that we can link to changes in per cent anisotropy. Unlike Ngatamariki (with near 100 per cent reinjection (Clearwater et al. 2015)), approximately 75 per cent of the geothermal fluid produced at Rotokawa is reinjected (Hernandez et al. 2015).
Per cent anisotropy at Rotokawa remains slightly more constant throughout the study period than the changes observed at Ngatamariki. No changes can be readily associated with changes in injection (Fig. 10). A slight downward trend of k cum in the later part of the catalogue appears to be negatively correlated with a much stronger increase in v P /v S ratios. As discussed in Section 1.2, normally a correlation between v P /v S and per cent anisotropy would be expected for liquid filled cracks (Unglert et al. 2011). This suggests that a high proportion of cracks within the reservoir are unsaturated. The fluid within the reservoir is known to be two phase so steam filled unsaturated cracks are present (McNamara et al. 2016). The majority of the variation in v P /v S ratios observed in Rotokawa is from P-wave velocities (we observe a higher correlation between v P and v P /v S than between v S and v P /v S ). This, along with the negative correlation between v P /v S and per cent anisotropy, is consistent with either unsaturated (dry or vapour-filled) cracks closing or unsaturated cracks becoming saturated (O'Connell & Budiansky 1974).
By expanding the catalogue to include the commissioning of the Nga Awa Purua power station in 2010 significant changes in production and injection would be captured. This would provide valuable information about the before state and how the anisotropy and v P /v S ratios react to the changes.

Pore-fluid pressure
We propose a threshold mechanism whereby if the pore fluid pressure drops below a certain value, cracks begin to close. Net volume is being removed from the south of Ngatamariki throughout the entire study period, however anisotropy either remains stable or increases until the increase in production (and proportionally smaller increase in injection) in late 2014 appears to trigger the rapid decrease in per cent anisotropy (Fig. 8). A drop in pore fluid pressure would also promote boiling causing the proportion of gas filled cracks to increase and thus decrease v P /v S . Boiling increases the pressure and so acts to stabilize the pressure drops, allowing the reservoir to reach a new equilibrium with a higher rate of production and a greater proportion of gas filled cracks in the reservoir.
While the cracks are aligned with S Hmax , the magnitude of S Hmin and pore fluid pressure (p f ) determine to which degree they are open or closed. The stress aligned cracks would begin to close as p f becomes less than S Hmin (Tod 2002). In northern Ngatamariki (with only injection) p f is greater than S H min , thus the cracks remain open. In southern Ngatamariki, if it is true that we are closing cracks and reducing pore fluid pressure, then we should also see a decrease in the rate of earthquakes despite an increase in production and injection in late 2014. This is exactly what occurs. Fig. 11 shows the cumulative number of earthquakes in southern Ngatamariki overlain by production and injection. There is a clear drop in number of earthquakes occurring beginning with the increase in production and injection in late 2014. The character of this seismicity is investigated by Hopp et al. (2019).
As we are in a predominately normal faulting regime we know that the maximum stress is vertical. Thus if a relationship between injection/production and pore-fluid pressure could be modelled then the absolute values of the principle stresses could be estimated based on the fact that we know when p f becomes less than S Hmin .

Implications for field development and operations
If the cracks we observe with shear wave splitting are related to those that facilitate fluid flow through the geothermal reservoir then results we find here may have implications in how geothermal reservoirs are managed. Our observations themselves demonstrate the value of seismic monitoring of geothermal fields before they are developed.
Firstly, it may be important to maintain high enough pressure in the reservoir such that σ /p f < 1 so that cracks remain open, allowing the geothermal fluid to flow. In a closed system such as Ngatamariki (Clearwater et al. 2015) this may be possible, however for systems where mass is lost (such as Rotokawa) it may still be achievable locally in areas where high permeability is important. Also, it may be important to maintain pore fluid pressure so that σ /p f is close to unity so that the rate of induced earthquakes is low. Monitoring of anisotropy with shear wave splitting could aid measuring and maintaining in this relationship.
Second, in finding the optimum position for new wells our results suggest that (at least in a normal faulting regime) to have a strong fluid connection wells should be placed along the axis of maximum horizontal stress. This is because these are the cracks that are easiest to keep open and appear to be the ones that respond most readily to changes in injection and production. Or, for when low fluid connectivity is required, wells should be placed along the minimum horizontal stress axis. This of course is with respect to local geology. Shear wave splitting fast polarizations may be useful in this situation for determining the local orientation of maximum/minimum horizontal stress to inform the placement of new wells, similar to what was previously suggested in other studies (e.g. Evans et al. 1995).

C O N C L U S I O N S
Shear wave splitting fast polarizations remain relatively uniform throughout the four year study period for both Rotokawa and Ngatamariki geothermal fields. Within Ngatamariki we observe mean fast polarizations that remain compatible with the NE-SW regional orientation of maximum horizontal stress and the strike of nearby faults. Within Rotokawa we observe rotations of the fast polarizations away from NE-SW towards N-S. While changes of this extent are not observed in the local estimates of horizontal stress orientation, they do support a local stress orientation within Rotokawa of less than 45 • (NE-SW) and show that strong local variation in the stress orientation is expected. We also observe correspondence between fast polarizations and the strike of the major Rotokawa faults, particularly in the southern portion of the field.
Shear wave splitting delay times and their associated per cent anisotropy show a strong relationship with production and injection within Ngatamariki. The beginning of field operations is marked by a higher than average per cent anisotropy. In the north of the field, where primarily injection occurs, increasing the injected volume corresponds to higher per cent anisotropy. In the southern and central portion of the field, where both injection and to a larger extent production occur, an increase in the total volume of fluid being removed is linked to below average per cent anisotropy, v P /v S ratios, and rate of earthquake occurrence. This observation is compatible with saturated cracks closing in response to a drop in pore fluid pressure. With no major changes in production/injection in Downloaded from https://academic.oup.com/gji/article-abstract/220/1/1/5561446 by Lawrence Berkeley Lab user on 03 March 2020 Rotokawa for our study period, changes in per cent anisotropy and v P /v S are consistent with either unsaturated cracks closing or unsaturated cracks becoming saturated.
Delay time tomography could be applied in the future to determine the spatial variation in percentage anisotropy similar to Johnson et al. (2011) or with a Bayesian approach. Expansion of the earthquake catalogue particularly to capture major changes in production/injection in Rotokawa would be useful in further understanding the relationship between anisotropy, cracks, and field operations (particularly in a field without 100 per cent re-injection). The modelling of pore-fluid pressure with respect to injection and production could also give an indication of the stress magnitudes in the reservoir (by comparison with crack closing/opening signals in per cent anisotropy).

A C K N O W L E D G E M E N T S
We would like to thank Mercury and the Rotokawa Joint Venture (Mercury and Tauhara North No.2 Trust) for their support of this research in awarding an MSc scholarship to SM, providing access to data, funding and granting permission to publish this work. We are also grateful to Steven Sherburn for providing the earthquake data and Emma Greenback for providing manually picked events from 2013. Stations WPRZ and ARAZ are operated by GeoNet, funded by the Earthquake Commission. We thank editor Donna Blackman and two anonymous reviewers and an anonymous associate editor for comments that helped to improve the paper.
Data processing and analysis were carried out by SM. Matched filtered events were provided by CH. SM and MKS performed the writing and editing with significant contributions from SMS and CH.
Figures and maps were made using the ggplot2 (Wickham 2009), ggmap (Kahle & Wickham 2013) and ggsn (Santos Baquero 2017) packages. All packages are within the R software environment (R Core Team 2013). Maps are sourced from the LINZ Data Service and licenced for reuse under the CC BY 4.0 licence. The MFASTR R package can be found on the github repository at git.io/vNXSX and contains the codes for conducting automatic shear wave splitting. MFASTR also includes the functions used for analysis and processing in this paper along with the summary files containing the raw shear wave splitting and traveltime data. Geosciences 2018, Napier, New Zealand, p.

S U P P O RT I N G I N F O R M AT I O N
Supplementary data are available at GJ I online. Figure S1. Station completeness for Ngatamariki stations. Each point represents a day where at least one class0 S-wave pick was made. The vertical blue line shows the beginning of production and injection. Figure S2. Station completeness for Rotokawa stations. Each point represents a day where at least one class0 S-wave pick was made. Figure S3. Difference between the three S-wave pick times in relation to the automatic P-wave pick. The AIC pick (red) is different between the two components. The pick time that we use for later processing is the earliest AIC pick (in this case on the north component). Figure S4. v P /v S for automatic versus manually picked events. Events are randomly selected from either 2014 (blue triangles) or 2015 (orange circles). Figure S5. Southern Ngatamariki (subset −38.58 < latitude < −38.55 and 176.18 < longitude < 176.20): subset of 577 events taken from the southern Ngatamariki cluster. Weekly mean per cent anisotropy (k, black/solid) and weekly mean v P /v S ratios (v P /v S , red/dashed) are plotted in both their cumulative (acting as a lowpass filter) and unfiltered forms. Daily injection/production (tonnes) is plotted in blue/dot-dashed and green/dashed (with a 7-d median filter applied). 95% confidence intervals are indicated by dotted lines and the vertical blue/dashed line marks the beginning of production. Labels on the x-axis mark the start of that year. Figure S6. Northern Ngatamariki (subset −38.55 < latitude < −38.53 and 176.17 < longitude < 176.19): subset of 217 events taken from the northern Ngatamariki cluster. Weekly mean per cent anisotropy (k, black/solid) and weekly mean v P /v S ratios (v P /v S , red/dashed) are plotted in both their cumulative (acting as a lowpass filter) and unfiltered forms. Daily injection (tonnes) is plotted in blue/dot-dashed (with a 7-d median filter applied). 95% confidence intervals are indicated by dotted lines and the vertical blue/dashed line marks the beginning of production. Labels on the x-axis mark the start of that year. Table S1. Table of spickerC parameters used in that have changed from Castellazzi et al. (2015, supp. info.). Includes descriptions after Castellazzi et al. (2015, supp. info.). Table S2. Filters tested in MFAST processing along with the per cent of high-grade measurements that were made with that filter. Table S3. Summary of mean shear wave splitting parameters for each station from all high-grade events weighted by final grade. Includes standard errors (S.E.) and counts of events for each grade (#). Includes stations that pass a weighted Rayleigh test (p-value < 0.05) and only events with an angle of incidence less than 35 • . Table S4. Summary of mean shear wave splitting parameters for each station from all high-grade events weighted by final grade. Includes standard errors (S.E.) and counts of events for each grade (#). Includes stations that fail a weighted Rayleigh test (p-value ≥ 0.05) and thus does not have a meaningful unimodal mean. Only includes events with an angle of incidence less than 35 • . Table S5. Estimated mean weighted per cent anisotropy (k) and differential shear wave anisotropy (SW A) for events measured at each station.
Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.