Relocating microseismicity from downhole monitoring of the Decatur CCS site using a modified double-difference algorithm


 The injection of CO2 at the Decatur carbon capture and storage site has generated significant microseismic activity, which occurs in distinct spatial clusters up to approximately 2.2 km from the primary injection well. Accurate and precise event locations are vital for the characterization of the microseismicity to help understand the reservoir response to the CO2 injection, whilst enabling the identification of minor faults and fractures below the resolution of conventional active seismic imaging. However, microseismic monitoring of fluid injection sites, such as Decatur, is often performed using a network of borehole sensors often from a single well. While these downhole sensors have excellent detection capabilities, their poor azimuthal coverage limits the ability to precisely determine event locations. We have developed a modified double-difference relocation algorithm suitable for both 1-D and 3-D velocity models, and which incorporates differential back azimuth observations to allow the benefits of the original double-difference algorithm to be applicable to a downhole microseismic monitoring setting. Applying the modified double-difference algorithm to the microseismicity at Decatur, we have successfully relocated 4293 events. The relocation included over 59 million observations for 757 285 event pairs, split across seven geographic regions. Despite the majority of observations being recorded in only two boreholes, with an almost identical azimuthal coverage, the results have shown to be reliable with significantly reduced residuals and low uncertainties associated with the final locations. We have analysed the residuals in terms of their association with each geographic region, data type, station and individual events, to fully appreciate their influence in the inversion and the fit of the data to the final set of event locations. For each region, the relocated seismicity has become less diffuse with improved clustering, and with newly visible linear features often orientated in a NE–SW direction. These results show the potential improvements that can be made to microseismic event locations recorded by a borehole network with a limited and variable azimuthal distribution.


I N T RO D U C T I O N
The Illinois Basin-Decatur Project (IBDP) is the first demonstration scale project from the U.S. Department of Energy's Regional Carbon Sequestration Partnership program (Finley et al. 2013;Finley 2014). This program, run by the Midwest Geological Sequestration Consortium (MGSC) aims to explore the technical and economic feasibility of using geological formations for long-term CO 2 storage. During a first phase of injection, one million tonnes of supercritical CO 2 were injected between November 2011 and November 2014, at a depth interval of 1.92-1.94 km below sea level (BSL), within the Lower Mount Simon Sandstone. In April 2017, as part of the Illinois Industrial Sources Carbon Capture and Storage (IL-ICCS) project, a second phase of injection began 1.1 km to the northeast of the original injection well. This injection, which continues as of 2020, occurs at the depth interval of 1.81-1.87 km BSL in another high permeability and porosity section of the Lower Mount Simon.
The intracratonic Illinois Basin spans much of Illinois, also extending into parts of Indiana and Kentucky, and is dominated by a strike-slip stress regime with the maximum horizontal stress orientated approximately northeast southwest (Lahann et al. 2017). While natural seismicity is associated with the southern, and to a much lesser extent the northern part of the state, the central Illinois region, which incorporates the Decatur CCS site, is rarely associated with felt seismicity . The Precambrian basement at the Decatur CCS site comprises a highly fractured rhyolite. Above the basement, the Argenta formation, is characterized by a low porosity (< 10 per cent) and low permeability < 36 mD (Freiburg et al. 2016) sandstone. Above it is the Mount Simon sandstone reservoir itself, which is relatively heterogeneous. In the Lower Mount Simon, where both CO 2 injections are located, the porosity is highest, with an average of 22 per cent and 28 per cent for the CCS1 and CCS2 injection intervals, respectively. For the permeability, this averages 200 mD (with a maximum of 1000 mD) and 250 mD (with a maximum of 1066 mD) for the respective injection intervals (Freiburg et al. 2014;Williams-Stroud et al. 2020). Throughout the Mount Simon there are discontinuous lenses of low permeability mudstone baffles 1.5-3 m thick. Such a mudstone baffle occurs directly above the CCS1 injection interval, limiting the upward vertical pressure during the first phase of injection . Above the reservoir, the Eau Claire Formation is characterized by an extremely low permeability shale near its base, which acts as an effective cap-rock for the reservoir formation.
Microseismic monitoring at the site has primarily been performed using a network of 3-and 4-component geophones deployed in multiple boreholes, but has also been supplemented by a network of surface seismometers and accelerometers. Between December 2011 and July 2018, over 19 000 detections were made on the borehole sensors (Bauer et al. 2019). 5 397 of these detections were considered to be of sufficient quality to locate.
Analysis of microseismic activity at CO 2 storage sites is a key tool for understanding how the subsurface is responding to the injection of CO 2 . Induced seismicity is an indispensable imaging resource for long-term CO 2 storage and for monitoring its conformance, with the potential to identify the generation of pressure pathways caused by fluid injection operations. This is vital for maintaining the integrity of the storage site, optimizing and managing the injection programme, and mitigating seismic hazard.
While microseismic event locations have the potential to track the propagating pressure front caused by CO 2 injection, observed seismicity often reflects a more complex pattern of changing stress conditions. Pressure changes can occur far beyond the extent of the CO 2 plume, while stress and strain changes can occur at even greater distances from the volume of pressure change (Rutqvist 2012). This may result in seemingly disparate event clusters from which an interpretation must be made.
Accurate and precise event locations provide the basis for further interpretation of seismicity, and while event characterization (e.g. source mechanism determination, b-value analysis, and stressdrop) can provide valuable information for aiding geomechanical and fracture modelling (e.g. Bachmann et al. 2012;Goertz-Allmann & Wiemer 2013), reliable locations remain a prerequisite. By identifying fractures and faults from microseismic event locations, we are better able to identify cap-rock failures (or risk of failure) and potential pathways for leakage of CO 2 . The spatial distribution of these fractures helps us to understand fluid migration and stress transfer within or around the reservoir, providing insight into the behaviour of the injected CO 2 , the reservoir and the seismic hazard. Moreover, although more advanced active seismic processing utilizing diffractions can identify structural features at a resolution below the seismic wavelength (Khaidukov et al. 2004;Moser & Howard 2008;Schwarz & Krawczyk 2020), microseismicity may similarly reveal fractures that can not be identified from conventional active seismic imaging, providing further insight into the CO 2 storage site.
Determining accurate and precise event locations, however, is not straightforward, and is reflected in the multitude of location methods that exist. These can include both traditional arrival-time picking approaches (e.g. for an overview see Wuestefeld et al. 2018)  and waveform based approaches (e.g. Gharti et al. 2010;Chambers et al. 2014;Li et al. 2020). Event location using arrival-time picks is a versatile method applicable for numerous sensor networks, but requires precisely picked phases. In contrast, waveform based methods remove the need for phase picks and can perform well for events with low signal-to-noise ratio (SNR) by exploiting waveform coherency across a network. However, optimal performance requires a dense network of sensors and in the case of borehole networks, where waveform complexity is more apparent, secondary phases or coherent noise can generate spurious event locations. While both arrival-time picking and waveform based methods typically operate on single events to determine absolute event locations, an additional class of location algorithm focuses on resolving relative hypocentres for multiple events. These include the master event method (e.g. Ito 1985;Zollo et al. 1995), the double-difference (Waldhauser & Ellsworth 2000) algorithm, and the source-specific station term method (Richards-Dinger & Shearer 2000), which all aim to reduce systematic biases in arrivaltimes resulting from 3-D velocity variations that are unaccounted for. In applying such a relative location algorithm, the effect is often an improvement in the pattern of seismicity, where underlying structural features are revealed, due to the improved location precision.
Regardless of whether an absolute or a relative location procedure is applied, having a good spatial distribution of sensors greatly benefits the event locations. Good azimuthal coverage can be achieved using a network of surface sensors, but at the expense of detectability and depth resolution. Since we aim to understand the reservoir response to CO 2 injection at the smallest possible scale, to build the most complete picture of the subsurface processes, the detection of low magnitude microseismic events is required. For example, to observe displacements on the mm-to-m scale, and subsurface rupture lengths of metres to tens of metres, requires the detection of events with magnitudes between -1 and 1 (e.g. see Mc-Garr & Fletcher 2003;Boettcher et al. 2009;Kwiatek et al. 2011;Göbel 2014;Thingbaijam et al. 2017;Williams-Stroud et al. 2020). Improved detectability is often achieved by deploying a string of sensors within deep boreholes close to the injection interval. In addition, the vertical span of sensors aids to constrain the depths of event locations. However, although the particle motion from such borehole sensors may be used to constrain event azimuths (e.g. Oye & Roth 2003;Jones et al. 2014), without the availability of sensors from multiple boreholes, the azimuthal coverage will be insufficient for precise event locations. While including additional phase arrival information may significantly reduce the overall event location uncertainty (e.g. Goertz-Allmann et al. 2017b), accurately identifying the correct phases and modelling their traveltimes can be challenging, depending on the subsurface geology (e.g. Dando et al. 2019). Moreover, the injection of CO 2 will alter the velocity structure (Wuestefeld & Weinzierl 2020), which has the potential to affect event location accuracy depending on the volume of CO 2 , and its placement relative to the source-receiver paths. Ultimately, the aim is to derive event locations that are less sensitive to spatial and temporal velocity variations, and where the precision matches that of relative location methods, despite constraints imposed by a sensor network with limited azimuthal coverage.
In this paper, we demonstrate how to improve microseismic event locations recorded as part of the IBDP and the IL-ICCS project, using a modified double-difference algorithm. We show that by adding differential back azimuth observations into the inversion, we are able to determine event locations with improved precision, which are stable even in the case of poor azimuthal sensor distribution. We carefully analyse the resulting traveltime and back azimuth residuals to demonstrate the influence of different observation types, stations, events and geographic regions on the inversion, highlighting the care required when running such a large-scale inverse problem.

CO I N J E C T I O N A N D M I C RO S E I S M I C M O N I T O R I N G AT T H E D E C AT U R C C S S I T E
Microseismic monitoring at Decatur consists of sensors deployed in multiple boreholes, operational one and a half years before the start of the first injection (Fig. 1). Two four-component geophones were installed in the first injection well (CCS1) at depths of 1660 and 1493 m (named PS3-1 and PS3-2, respectively), while a separate geophysical monitoring well (GM1) 56 m to the northwest of the injection well was drilled and installed with 31 three-component geophones-29 sensors deployed between 418 and 844 m depth, and two additional sensors close to the surface at 97 and 164 m above sea level (ASL); the approximate elevation at Decatur being 206 m ASL. The sensors from these two almost colocated wells (CCS1 and GM1) were operational throughout the monitoring period and were supplemented by 5 three-component geophones that were first deployed in a verification well (VW2), 1.7 km north of the first injection well at depths of between 1659 and 1903 m BSL, before being moved to a second geophysical monitoring well (GM2) at depths of between 575 and 821 m BSL, located 1.1 km to the northeast of the first injection well. Monitoring from the VW2 and GM2 wells was inactive for much of the recording period, substantially reducing the azimuthal coverage.
For the period between December 2011 and July 2018, 5 397 microseismic events (Fig. 1) were located by applying a modified Geiger location algorithm Geiger (1912) that uses both P-and Swave first arrival times recorded on the available borehole sensors and polarization analysis to constrain the event back azimuth . The polarization analysis was limited to the PS3 sensor observations and if available, the VW2 sensors (Schlumberger 2019, personal communication). 4848 of these located events occurred during the first injection phase, with the remaining 549 events occurring in both the post-CCS1 injection period (from December 2014 onwards) and after the commencement of the second injection phase in CCS2 (from April 2017). Events occur in distinct spatial clusters with a lateral distance of up to 2.2 km from the CCS1   Table 1. For each box plot, the boxes represent the interquartile range (IQR) for a given month, the blue line within each box represents the median value, while the whiskers represent the upper and lower data limits within 1.5×IQR of the upper and lower quartiles, respectively. Data beyond the whiskers are considered outliers and are plotted as individual blue dots. injection well. Focal depths are distributed from the lower reservoir formations (Lower Mount Simon and Argenta) into the upper Precambrian basement. However, measurement errors of seismic arrival times and velocity model inaccuracy, especially from applying a 1-D velocity model, may lead to significant location uncertainties for individual events. Unfortunately, neither event location uncertainties, nor the arrival time data required to compute them, were provided in the original event catalogue. Table 1. Summary of the operational status of sensors from November 2011 to July 2018. The total number of sensors was 43, but the maximum operational at any one time varied between 26 and 38. All available sensors are shown in Fig. 1. Small variations in data continuity persist throughout the recording period, thus the exact station continuity is not represented in this  A timeline showing the distribution of events is shown in Fig. 2. The magnitudes range from -2.1 to 1.2, based on the catalogue data. The complete distribution is shown in Fig. 3.
In addition to the borehole stations shown in Fig. 1, the network was further expanded to include surface and near-surface stations. From 2013, the United States Geological Survey (USGS) installed a total of 7 shallow borehole stations comprising a 3-component (3C) geophone at 150 m depth supplemented by a 3C surface accelerometer. In addition, the USGS installed 10 surface stations, each combining a 3C broad-band seismometer and 3C accelerometer. In 2013, the Illinois State Geological Survey (ISGS) installed a further 5 3C seismometers in surface vaults. The locations of these surface and near-surface stations are shown in Langet et al. (2020), but were not used to locate the microseismicity shown in Fig. 1. From July 2013 to February 2015, the USGS located 179 events (Kaven et al. 2015) using their surface network. Note that only a small subset of events recorded on the downhole network are detected at the surface-approximately 9 per cent between July 2013  to February 2015. However, these surface sensors have proved particularly useful for determining focal mechanism solutions , and aided fault characterization and interpretation of reservoir processes . Combining borehole and available surface network data greatly improves the azimuthal coverage and therefore leads to better constrained absolute event locations . For 23 events located in the northernmost cluster, Langet et al. (2020)    Number of events and corresponding number of event pairs and input observations that were used for the relocation in each region. These numbers reflect the final relocated solution of each region meaning observations and events, which have been removed between iterations (e.g. due to weighting constraints), are not included. Two events were shared between region A and region B. Only the region B solutions for these two events are shown in the subsequent results, but they are included in the numbers given for region A in this and that the 23 events are of sufficient quality to determine focal mechanism solutions, these uncertainties can be considered the best-case estimates for events within this cluster. For the catalogue data, based solely on the downhole observations with highly variable signal-to-noise ratios and waveform characteristics, location uncertainties are expected to be much larger.

A P P L I C A B I L I T Y O F T H E D O U B L E -D I F F E R E N C E A L G O R I T H M T O D O W N H O L E M O N I T O R I N G O F M I C RO S E I S M I C I T Y
The relative relocation of seismicity using the double-difference algorithm of Waldhauser & Ellsworth (2000) is well established for improving the relative locations between events. The methodology uses P-and S-wave differential traveltimes for pairs of events recorded at common stations. These observations can be derived directly from phase picks or from time delays based on waveform cross-correlation. The methodology assumes that the ray paths to a common station will be similar for small interevent distances compared to the event-station distance. Further, it is also assumed that any velocity heterogeneity along the ray paths will be larger than the interevent distances, although this requirement depends on the location of the heterogeneity along the ray path (i.e. smaller heterogeneities closer to the station will not change the relative ray paths between the two events). With similar ray paths for an event pair, any traveltime differences can therefore be largely attributed to changes in the spatial offset between the events. An overview of the original double-difference method is given in the Appendix.
The application of this double-difference relocation methodology is often focused on local and regional seismicity (e.g. Waldhauser & Ellsworth 2000;Yang et al. 2005;Chiaraluce et al. 2011;Sippl et al. 2013), including cases of induced seismicity (e.g. Albaric et al. 2014;Woo et al. 2020), with teleseismic examples for the relocation of global catalogues also existing (Waldhauser & Schaff 2007;Pesicek et al. 2010;Diehl et al. 2013). These studies are characterized by having a good azimuthal distribution of stations, critical for well-constrained event locations. Microseismicity triggered or induced by fluid injection, for example, at sites associated with carbon capture and storage (CCS), hydraulic fracturing, or enhanced geothermal systems, is less likely to be located using a seismic network with optimal geometry. Instead, these sites often have a monitoring network comprising one or more borehole strings, instrumented with 3-component geophones. For such monitoring geometries, the microseismic locations are often constrained by using the polarization angles derived from the 3-component data, in addition to the traveltimes (e.g. Oye & Roth 2003;Jones et al. 2014). Such polarization data are not included in the double-difference algorithm. Bai et al. (2006) analyse the effect of different station distributions on the location error for the double-difference method. They show how network geometries with high azimuthal gaps produce unstable double-difference relocations. To achieve a 1 km relative hypocentral error or a 0.5 km relative epicentral error to the 95 per cent confidence level, Bai et al. (2006) recommend that at least 15 ± 2 stations are installed within 100 km of the epicentre, and that the azimuthal gap is less than 210 • ± 15 • -conditions that are clearly not met with a single borehole string. In the case of multiple, well distributed borehole stations, the application of the double-difference method will, however, still be viable, as demonstrated by Castellanos & van der Baan 2013).

M O D I F I C AT I O N O F T H E D O U B L E -D I F F E R E N C E A L G O R I T H M
To account for the lack of the azimuthal constraint at the Decatur CCS site, where only the GM1 and PS3 stations are continuously operational, there is a clear need to modify the original doubledifference algorithm for it to generate stable relative locations. Previous work to address this issue has involved a variety of approaches. Li et al. (2013) presented a methodology that solves for both location and anisotropic velocity structure by modifying the double-difference tomography method of Zhang & Thurber (2003 in a synthetic study. The modification includes incorporating differential back azimuth observations into the system of linear equations providing additional constraints to the lateral perturbations between event locations. The method has subsequently been demonstrated on a small data set of 68 microseismic events (Li et al. 2014). In contrast, Tian et al. (2016) reformulated the relocation problem in 2-D to eliminate the need for an azimuthal constraint, by solving simply for the change in horizontal distance, in addition to depth and origin time. The relocated events are then projected back into 3-D space by using individual event back azimuth observations from P-wave polarization analysis. The 2-D formulation only requires a simple 1-D velocity model and provides a significant computational benefit. While 1-D velocity models are much more frequently used in microseismic monitoring studies than 3-D models, the method is inapplicable if a well-constrained 3-D velocity model is available. Furthermore, unlike in the method of Li et al. (2013), where the differential back azimuth observations reduce the relative azimuthal uncertainty between event pairs, the single  Zhang & Thurber (2003), to simultaneously solve for both locations and velocity structure, with application to downhole microseismic monitoring. To provide an azimuthal constraint to the locations, they use the polarization data to obtain back azimuth observations for each event. These single event observations are then used as a constraint in the inversion. However, similar to Tian et al. (2016), these observations are still incorporating the errors associated with individual back azimuth measurements.
We based our approach to the relocation of the microseismicity at the Decatur CCS site on the methodology by Li et al. (2013), but limited the inversion to the relocation problem, avoiding the additional complexity associated with simultaneously resolving for changes in the velocity structure. Thus, the differential back azimuth residuals can be defined for an event pair (i, j) as where φ obs and φ theo are the observed and theoretical back azimuth, respectively, at station k. Similar to eq. (A7), the differential back azimuth residuals can be linearly related to the hypocentral perturbations for an event pair, using For a 1-D velocity model, the partial derivatives for the back azimuth observations can simply be expressed as Incorporating the back azimuth observations into the data vector in eq. (A9), we can expand our inversion scheme to Here, d t and d φ are the data vectors containing the differential traveltime residuals and differential back azimuth residuals, respectively. Similarly, the partial derivatives have been split into matrices containing those of the traveltime sensitivities, G t , and those of the back azimuth sensitivities, G φ . Also of note is the diagonal weighting matrix W, which differs depending on the data type, that is W t or W φ . In fact, since the traveltime data itself can consist of both P-wave and S-wave derived traveltime data, either from picks (i.e. t Ppick or t Spick ) or from cross-correlation time-lags (i.e. t Pcc or t Scc ), the system of equations can be further divided such that ⎡ Furthermore, we apply three additional constraints following Waldhauser & Ellsworth (2000) that are user-controllable. First, since the double-difference method resolves relative locations, we can constrain the mean shift of all N events to be zero. Where we have confidence in the absolute location of the event cluster as a whole, as expressed via the starting locations, this constraint will aid the stability of the solution. The constraint is implemented by extending eq. (8) by an additional four equations such that Secondly, since the G matrix is extremely sparse with only 8 nonzero elements in each row, corresponding to the event pair for a given observation, we apply a scaling by normalizing the L2-norm of each column to improve the numerical stability. Finally, we apply a damping factor to the inversion to further reduce the likelihood of instability in the case that events are poorly linked. However, such regularization is regarded as secondary to careful data selection as set out in the subsequent section. In selecting a damping factor, multiple trial runs were completed, to ensure optimal data fitting to the model parameters, whilst retaining stability. Given the sparsity of the G matrix and the number of events and observations recorded at the Decatur CCS site, we solve for the model parameters in eq. (8) using the LSQR algorithm (Paige & Saunders 1982), iterating the inversion until the rms residual reduction is negligible.
For the different weighting factors shown in eq. (8), we implement a weighting scheme similar to Waldhauser & Ellsworth (2000). A priori weights are defined for each station, data type, phase and iteration. However, to ensure that the residuals derived from arrival-time picks, cross-correlation time-lags and back azimuths are balanced with respect to each other, we normalize the a priori weights associated to each of these data types using the rms of the data vector for that particular data type. For the cross-correlation derived data, we multiply these a priori weights with the squared cross-correlation coefficient. The inversion is run for multiple iterations with this weighting scheme until a stable solution is produced. Subsequent iterations are then subjected to an optional reweighting based on the residual misfit, the interevent distance or both. For this study, no reweighting of the back azimuth observations was applied, with only the a priori weights implemented for this data type. For the arrival time measurements and cross-correlation data, reweighting by residual misfit provided the most stable solutions.
While we have defined the partial derivatives for a 1-D velocity model in eqs (3)-(6) and (A2)-(A5), which is applicable for the Decatur CCS site, our implementation also includes the definitions required for implementing a 3-D velocity model. Both sets of definitions are dependent only on station-based traveltime tables, which can be computed prior to the inversion, avoiding costly ray tracing at each iteration.

Arrival-times, back azimuths and quality control
Despite obtaining the data, locations, origin times, and magnitudes associated with the 5397 microseismic events shown in Fig. 1 (Bauer 2019, personal communication), the arrival-time and back azimuth data used for deriving these locations were unavailable. Careful data selection is critical for reliable solutions to inverse problems. In consideration of variable waveform characteristics (Goertz-Allmann et al. 2017a) and inconsistent data quality, the complete data set was therefore manually picked to ensure reliable and high quality data as input for the double-difference relocation. Out of the 5397 events previously located, we were able  to reliably identify phases and measure arrival times for 4320 events. Only seven sensors had temporally consistent sensor orientation information, such that they could be used to derive back azimuth observations. These correspond to the two PS3 and the five VW2 sensors. Prior to picking, the data were first filtered with a 60 Hz notch filter with a bandwidth of 15 Hz to remove electrical noise on many sensors, followed by a 40-200 Hz bandpass filter.
An overview of the number of P-wave picks, S-wave picks and back azimuth observations for the 4320 events is shown in Fig. 4. The red lines show the maximum number of possible observations for a given time period, which varies due to the operational status of different stations. A summary of when stations were active or inactive is shown in Table 1. From Figs 4(a) and (b), it is recognizable that the number of S-wave picks per event is consistently higher and less variable than the number of P-wave picks per event. This is due to the higher signal-to-noise ratio of the S-wave. In total, 114 530 P-wave arrivals were picked compared to 124 446 S-wave arrivals, with a median of 27 and 30 picks per event, respectively. There exists a notable change after the end of the CCS1 injection (Figs 4a and b), both in the number of picks per event and to a less extent in variability (e.g. June 2015 and January 2016). Whilst there are much fewer events in the post-CCS1 injection period, rendering these observations less statistically significant, the decrease in pick count reflects a much lower signal quality that was observed for these events. For the back azimuth observations shown in Fig. 4(c), we consistently obtained two observations per event (i.e. from both of the PS3 sensors) until August 2013, after which the VW2 sensors become active, providing additional measurements. Post-CCS1 injection (i.e. after November 2014), there is a reduction in the number of observations per event, which is followed by VW2 becoming inactive in March 2015. Thereafter, we consistently achieved only one back azimuth observation per event, due to much poorer data quality on one of the PS3 sensors. Despite this, we were still able to achieve stable inversion results for these events. In total, 13 902 back azimuth observations were made.
To perform the double-difference relocation, the 4320 events were divided into seven regions labelled A to G, based on their spatial distribution, which are assumed to be structurally independent (Fig. 5). The events within each region are then relocated separately. Relocation using multiple regions was critical for two main reasons. First, the relationship between individual event pairs and between event pairs and receivers (e.g. number of observations, interevent distances, station-to-event distances, etc.) varies significantly geographically due to the large area over which microseismicity is occurring. This requires different data selection criteria and regularization of the inversion to ensure a stable and reliable solution. Secondly, the large number of microseismic events provides an additional restriction in terms of computational demand, and the simultaneous relocation of all events would prove computationally challenging without using high performance computing facilities.
To ensure the input data to the inversion is of a high quality, outliers in arrival time picks and back azimuth observations must be removed whenever possible. This is achieved in several steps. First, all traces were visually inspected as part of the manual picking process, and phases were picked as consistently as possible. When both head waves and direct waves were present, the first arrival was picked. This is consistent with the theoretical traveltimes we computed. For the back azimuth observations at the PS3 and VW2 sensors, the P-wave particle motion used to derive the back azimuth was visually inspected on the hodograms to ensure linearity. Secondly, individual event traveltime and back azimuth residuals based on the starting locations were analysed. Outliers were inspected and either repicked, removed, or confirmed to be valid observations. In the case of the more error-prone back azimuth observations, any remaining residual greater than ±20 • resulted in the removal of that observation. The resulting set of individual event residuals [( t obs − t theo ) and (φ obs − φ theo )] are summarized in Fig. 6. All remaining outliers represent either a misfit in the original locations (e.g. due to velocity heterogeneities not accounted for in the velocity model) or from mislocations introduced by errors in the original (and unavailable) observations. Across all events, the median P-wave residual is approximately 0.0 ms, the median S-wave residual is -4.7 ms, and the median back azimuth residual is -1.16 • . Despite the large reduction in number of events after the end of the CCS1 injection and a general degradation in data quality, the residuals in Fig. 6 show fewer outliers and slightly less variance for this time period. However, this can be attributed to an overall reduction in number of observations, as shown in Fig. 4. Splitting these residuals into their respective regions that we have defined for performing the doubledifference relocation (Fig. 7), we can see that there are no clear systematic differences between the event residuals for each region. We do however, see slight differences in the number and range of the statistical outliers for each region. Region D shows the largest range in residual outliers for all three observation types, which is due to the significantly higher number of events and observations for this region. A similar observation can be made for region G, which has the second highest event count. Regions A and C are the most distant from the GM1 and PS3 sensors, which make up the bulk of the observations. This is perhaps reflected in the higher S-wave residual range for region A, and the higher P-wave residual range for region C. However, it should be emphasized that it is only the statistical outliers that are showing relatively large values, while overall the residuals are remarkably consistent. This suggests that our arrival-time and back azimuth measurements are likely very similar to the observations used to derive the original locations.

Data constraints and cross-correlation analysis
While individual event residuals may help to identify spurious arrival-time and back azimuth measurements, the double-difference algorithm itself relies on event pair observations. We therefore apply additional constraints to these observations to ensure that the double-difference assumptions are valid, event pairs are well constrained, and that the size of the input data does not become prohibitively large. For all regions, we require at least 16 traveltime observations and one back azimuth observation per event pair. In addition, we limit the event pairs based on interevent distance to help ensure ray path similarity. This varies from between 120 m for region D, located relatively close to the sensors and exhibiting a high number of events, to 400 m for region A, located at a greater distance from the monitoring network. Furthermore, the differential traveltime residuals were analysed for each region by comparing them with an expected delay time between the events based on a constant velocity to identify outliers; none were found. For the differential back azimuth residuals, the same constraint as for the individual residuals of ±20 • was applied for each event pair observation.
To determine the differential traveltime residuals derived from cross-correlation time-lags, we use only the event pairs and the observations that remain after applying the aforementioned data constraints. Although observations from event pairs with high interevent distances will be uncorrelated and thus down-weighted during the inversion, we prefer to improve efficiency by limiting the input data, whilst being confident that we are only selecting data that will positively contribute to the relocations. The P-wave data are Downloaded from https://academic.oup.com/gji/article/227/2/1094/6315329 by guest on 05 August 2021 Figure 15. The percentage rms residual reduction for each of the five data types (columns) and for each of the seven regions (rows A-G). The absolute rms reduction values are given in Table 3. cross-correlated using a 40 ms time window starting 10 ms before the P-wave arrival, while the S-wave cross-correlation uses a 50 ms time window, starting 10 ms before the S-wave arrival. Each component of each sensor is cross-correlated separately. However, since the orientations for the GM1 and GM2 sensors are both unknown and temporally inconsistent, we additionally compute cross-correlations using the geometric mean of the horizontal components. A final set of time-lags and cross-correlation coefficients are subsequently compiled by selecting the values for each individual event pair with the highest cross-correlation coefficient across all components (including the geometric mean). A summary showing the distribution for each region of both the cross-correlation coefficients and the time-lags are provided in the supplementary material. For the Pwave data, the median cross-correlation coefficient varies between Table 3. Rms residual reduction for each data type and region. The percentage reduction is shown in parentheses. A visual representation is given in Fig. 15 a minimum of 0.67 in region D to a maximum of 0.81 for region G. The interquartile range (IQR) is also remarkably consistent, varying between 0.16 in region D to 0.19 in region B. For the S-wave data, we observe both higher median cross-correlation coefficients and smaller statistical dispersion for all regions. The median S-wave cross-correlation varies between 0.80 (region F) and 0.85 (region G), while the IQR varies between 0.08 (region F) and 0.11 (region C). The cross-correlation derived time-lags for both the P-wave and S-wave data are also relatively small and with low statistical dispersion for all regions, providing further confidence in the manually picked arrival-times. The median time-lag for each region does not exceed 1 ms for either the P-wave or S-wave data, while the IQR for the P-wave time-lags varies between 2 ms (region G) and 8 ms (region D), with the IQR for S-wave time-lags varying between 4 ms (region D) and 7 ms (region A).
In addition to ensuring that the observations from both the arrivaltime derived and cross-correlation derived differential traveltimes match with regards to identical event pairs and stations, we provide the constraint that all events must have both a back azimuth observation and traveltime measurements, such that all model parameters are accounted for. The final constraint involves the clustering approach implemented in the original double-difference (HypoDD) software (Waldhauser & Ellsworth 2000). This clustering is based on how well event pairs are interconnected based on the number of traveltime observations, which we achieve using an agglomerative hierarchical clustering algorithm (MATLAB 2019). Where events are poorly interconnected, they are split into separate clusters and relocated separately. Out of the seven regions, only the processing of events from region D resulted in further subdivision due to poor interconnectivity. Out of the three clusters of events generated for region D, two contained only two events and were subsequently excluded from the analysis. The final set of events that was relocated totaled 4293 of the original 4320 micro-earthquakes. A summary of the number of events and observations used in the relocation for each region, is shown in Table 2. A visual representation of the observations for each event pair is available in the supplementary material. Although the total event count for region A is 138 in Table 2, this includes two events that were also relocated as part of region B due to a deliberate overlap in both regions. For the final results we assigned these two events to region B, since there was a larger rms reduction of the events' residuals for this region. The greatest change in location occurs for region A with a median displacement of 151 m. For this region, located furthest away from the monitoring network, we observe several interesting features. First, in map view (Figs 8 and 11), both the spatial distribution and overall orientation of the seismicity has significantly changed. A larger central cluster has been generated with a smaller separate cluster identifiable to the southeast. While outliers still exist, the orientation of the clusters is now much closer to the critical stress direction in relation to the maximum horizontal stress direction. In the east-west cross-sections (Figs 10 and 13), the starting locations for region A show that the events are split into two separate groups-one within the Argenta, the other within the Precambrian basement. This is indicative of the effect of a layered velocity model, generating artificial layering within the individual event locations. For the relocated events, this effect has been removed, with almost all events being located within the Precambrian basement.

Double-difference relocation results
For region B, the total change in location is secondary to region A, with a median shift of 73 m. Although there is still appreciable scatter in the event locations, three separate clusters emerge, which were previously indistinguishable. The easternmost of these is clearly visible in the east-west cross-sections (Figs 10 and 13), where the lateral scatter has been minimized.
The northernmost region, C, exhibits lower overall displacements compared to all other regions, with a median shift of 37 m. However, in contrast to the other regions, depth is the highest contributing component to these location changes. In map view (Figs 8 and 11), the scatter to the northeast and east has reduced, and the seismicity is much more tightly clustered. In the north-south cross-sections Figure 16. The percentage rms residual reduction at each of the 43 stations (columns), for each of the five data types (subplots), and for each of the seven regions (rows labelled A-G). Stations with no observations for a particular data type and region are shown as black. The colour scale saturates at -100 per cent. (Figs 9 and 12), two separate clusters are recognizable-a larger, steeply dipping cluster, with a narrow subvertical linear feature within it, and a smaller, slightly deeper cluster, associated with lower magnitude events.
Region D is by far the largest region comprising 2299 events relocated based on over 41 million observations ( Table 2). The relocation displacements, however, are relatively small with a median of 38 m. This is perhaps due to its central location, with respect to the monitoring network, allowing for well-constrained initial locations. Nevertheless, events cluster more tightly after the relocation, best observed in Fig. 11. The westernmost cluster of region D (Fig. 8), in particular, is much more linear, and all clusters are now broadly aligned with the critical stress direction.
For region E with a median displacement of 67 m, we likewise observe an enhanced clustering of events, most apparent in the south of the region (Figs 8 and 11). Moreover, in the east-west crosssections (Figs 10 and 13), events are now more restricted in depth to both the base of Lower Mt. Simon and the Argenta formation, with fewer events above the injection depth or in the Precambrian basement.
Region F, located to the southeast and totalling 188 events, has a median relocation displacement of 58 m. Similar to region E, the initial locations are relatively diffuse, while the post-relocation distribution appears more clustered (Figs 8 and 11). In particular, the southernmost cluster appears more linear with an orientation closer to the critical stress direction in relation to the maximum horizontal stress direction. North of this cluster, the events still display some scatter, but there is a clear concentration of events that was not apparent prior to the relocation. In the vertical direction, there do not appear to be many significant changes, with a median change in depth of only 15 m.
Finally, for region G, which has the second highest number of events (903), we clearly observe a single event cluster both prior to and after the relocation (Figs 8 and 11). However, the event distribution becomes more well defined forming an 'S'-shaped structure, approximately 750 m in length. The structure itself is more detailed and features linear NE-SW substructural features. The hypocentres for this region extend from the Argenta formation into the Precambrian basement, where most of the seismicity occurs (Figs 12 and  13). In terms of changes in depth, there are minimal differences, with a median change of only 6 m (IQR of 9 m).

Residual analysis
To assess the quality of the inversion, the results need to be assessed in terms of both the residual reduction and the location uncertainties. While the location uncertainties will help visualize the individual uncertainties based on the data and method, analysis of the residual reduction will demonstrate the performance of the inversion in terms of fitting the model (i.e. the event relocations) to the data.
Since we are inverting five different sets of observational data, each with their own weighting, we obtain a better overview when we split the achieved rms residual reduction into their corresponding data types as shown in Table 3 and Fig. 15. Fig. 15 demonstrates that the cross-correlation derived S-wave observations consistently provide the greatest reduction in residuals, ranging from approximately 41 per cent for region D to 64 per cent for region A. This is partly influenced by the strong weighting of this data type and the high number of S-wave observations available. It is also notable that both the cross-correlation derived P-and S-wave arrival time residuals show a similar (although typically slightly higher) rms reduction as their non cross-correlation derived equivalents. This suggests that there is relatively good consistency between the data types, despite the benefits that the cross-correlation data provide. Finally, we observe that the back azimuth observations generate the Figure 19. The number of observations as used for the double-difference relocations for each event (columns), for each of the five data types (rows), and for each of the seven regions (subplots labelled A-G). smallest residual reduction. Since the number of observations for these data are only a fraction of the arrival-time data (Table 2) and are weighted only to ensure stability, this behaviour is as expected.
Whilst the breakdown of residuals into their respective data types and regions provides insight into the performance of the inversion, an analysis highlighting weaknesses in the data fit is beneficial to avoid overinterpretation and to identify potential data deficiencies, which were not previously evident. In Fig. 16, we show the residual reduction associated with observations for each station, divided into each data type and region. It is immediately clear that for all applicable data types, the observations from the five GM2 stations perform poorly, with significant increases in the final residuals relative to the starting residuals, particularly for regions D and G. The exception to this are observations for regions B and E, where GM2 observations show improvements in the residuals for the picked and cross-correlation derived S-waves as well as the cross-correlation derived P-wave (although not for GM2-3 for region B). Comparing the performance of each station in terms of residual reduction with Table 4. Variation in the rms of the scaled (but unweighted) final residuals that are used in the uncertainty analysis. The rms is shown for each input data type, and for each region.    the number of observations at each station (Fig. 17) reveals the disparity between GM2 and the other stations; for example, the number of observations derived from the GM2 sensors are many orders of magnitude below the GM1 station observations. These observations will therefore have a limited influence on the inversion. Regardless of this lack of influence, however, there still remains a fundamental conflict between the final relocations and the GM2 observations for regions D and G. In contrast, at GM1-30 and GM1-31, which have comparably few observations, there is generally a reduction (i.e. improvement) in residuals for the final relocations. Despite the problems associated with the GM2 data, it is reassuring that the VW2 observations, which are made to the north of the Decatur CCS site, compared with the dominant GM1 observations from the south, generally provide a good fit to the final solutions, with the exception of the picked P-wave for region A. Finally, the residuals can also be assessed on an event-by-event basis to identify particularly poorly performing individual events. This is shown in Fig. 18, with the corresponding number of observations displayed in Fig. 19. While there are events which show a net residual increase (i.e. deterioration) for certain data types (particularly picked P-waves, and to a lesser extent back azimuth observations), the cross-correlation derived observations always generate a positive residual reduction for both the P-wave and S-wave data, with the exception of two events in region A and four events in region D, where cross-correlation derived P-wave data show a deterioration. For the back azimuth data, while these data have been weighted relatively low, it is worth noting that there are still numerous events with deteriorating back azimuth residuals. For regions A-G, the proportion of events with negative back azimuth residuals total 12, 9, 16, 0.5, 0.8, 1 and 0.1 per cent, respectively. In total, across all data types and regions, there are 326 events (i.e. 7.6 per cent of all events) for which at least one residual data type deteriorates. Considering the general residual improvement for the cross-correlation data for these events, their removal would be excessive. Moreover, the elimination of these events does not alter the results to a degree that affects their interpretation (see Supporting Information).

Location uncertainties
To assess the uncertainties of the relocated events, we base our approach on the bootstrapping methodologies (Efron 1982) as outlined by Shearer (1997) and Waldhauser & Ellsworth (2000). Instead of applying a predefined distribution of errors (e.g. Billings et al. 1994;Wuestefeld et al. 2018) on the double-difference traveltime observations, we use the final residuals to provide an estimate of the observation errors. Subsequently, these residuals are randomly sampled with replacement (i.e. 'bootstrapped') and used as input data to the double-difference inversion to provide an estimate of changes in the final set of model parameters. This bootstrapping procedure is repeated 1000 times to generate a distribution of model parameter perturbations for each event from which location uncertainties can be derived.
The final residuals are computed by first predicting the data using the model parameters at the final event locations as well as the partial derivative matrix, The final residuals are then defined as the difference between these predicted data and the initial double-difference observations (d initial ), We apply an additional scaling factor (s) to account for the number of degrees of freedom (f) in the double-difference inversion using s = n/(n − f), where n is the number of residuals. For the traveltime differences based on either arrival-time picks or cross-correlation data, f will be eight, whilst for the back azimuth differences, f will be only four due to the lack of depth and origin time constraint for each event pair. With the large number of observations, even within each data type, this scaling factor will have a negligible effect (the maximum residual change is only 0.02 per cent, which is applied to the back azimuth residuals in Regions A and F), but has been included as a measure to ensure the methodology is applicable to smaller-scale problems in the future. In addition to the different scaling factors required for the various data types, the final residuals will have their own associated errors based on the data type, and should therefore be treated independently. For example, there will be inherent differences between the uncertainties based on the S-wave arrival times and the P-wave arrival times. Moreover, the back azimuth observation residuals can clearly not be substituted for traveltime residuals and vice versa. For these reasons, we treat each data type individually for the bootstrapping, such that the resampling only occurs within each data type. A summary of the residuals used for the uncertainty analysis is shown in Table 4. The choice of 1000 bootstraps for each data type is perhaps excessive, but repeated runs for region A using 500 bootstraps showed small variations in the location uncertainty for a single event, prompting the increase to 1000 bootstraps, which produced consistently stable uncertainties. Finally, 95 per cent confidence uncertainty ellipsoids were calculated from the eigendecomposition of the covariance matrix for the bootstrapped location results.
The final location uncertainties are shown in Figs 20-22. As expected, region A, which is most distant from the monitoring network,  Fig. 20 is the influence of the GM1 and PS3 observations, with a clear azimuthal dependence around these boreholes. Despite this influence, the overall orientation of the clusters in each region after the relocation appear valid, with the uncertainties having less of an effect than the maximum horizontal stress direction.

D I S C U S S I O N
The application of the modified double-difference method to the microseismicity at Decatur has demonstrated the validity of the method to a downhole monitoring setting, whilst providing valuable improvements to the event locations. Separating the data for the relocation into seven different regions has ensured stable and reliable solutions by tailoring the inversion parameters and regularization depending on the distribution of event pairs and their recording stations. Such an approach can be considered valid if the seismicity within each region occurs independently from the other regions. This is borne out in the relocation results, which show that the events within each region become better defined in discrete clusters, with very few events moving from one region to another.
Before commenting further on the relocation results, it is worth considering some additional challenges associated with applying the double-difference method to a microseismic monitoring setting, beyond the modifications that were required to aid the lack of azimuthal coverage. The velocity structure at a fluid injection site is often characterized by multiple layers representing the permeable reservoir in which the injection will take place, with impermeable layers above and below. Michelini & Lomax (2004) demonstrate how the use of an incorrect velocity model, for example assuming a gradient model in place of the more representative layered model (or vice versa), will lead to bias and errors in the relative locations. The assumed and modelled ray paths will be markedly different from the true ray paths due to the assumption, for example, of diving waves in place of interface waves. Such errors in the model will lead to incorrect partial derivatives in eq. (A7) and ultimately the errors and bias observed by Michelini & Lomax (2004). While there are still likely errors in the 1-D velocity model we have used (Fig. 23), and we do not explicitly use a layered model, which would best model any interface waves, there is still significant detail with minimal smoothing, which reduces the risk that ray paths are misrepresented.
A related but perhaps more often overlooked challenge is the potential breakdown of the assumption that closely located events will have similar ray paths, irrespective of the validity of the velocity model. In a layered or vertically heterogeneous velocity model, with subhorizontal ray paths to borehole sensors, small changes in location may result in very different ray paths (Dando et al. 2019). For the Decatur CCS site, this can be observed in Fig. 23, which shows significant differences in the ray paths for both the first arrival and the direct arrival for two separate event locations. Whilst these differences in the ray paths may result in changes to the signal characteristics that can be used to discriminate between event locations (Goertz-Allmann et al. 2017a), there is a breakdown in the ray path similarity assumption for the double-difference method. The consequence is that any observed relative traveltime residuals (eq. (A6)) for this event pair and either of the PS3 sensors in this example, cannot be assumed to be fully attributable to the focal region. However, in the case that the velocity model is well constrained and valid for the modelled ray paths, such a concern should not fundamentally limit the ability to accurately relocate the events. Moreover, having observations from multiple sensors for a large number of event pairs helps to accurately constrain the relocation results. It should also be noted that the event locations in Fig. 23 are deliberately chosen to demonstrate a potential breakdown of the ray path similarity assumption. For most parts of the Decatur field, ray paths of closely neighbouring events do not significantly change.
Since the injection of CO 2 will have an effect on the velocity structure (e.g. Wuestefeld & Weinzierl 2020) and the timing of the microseismicity occurs throughout two separate injection phases and a post-injection phase, there will naturally be unaccounted for changes in the 1-D velocity model used to locate the events in this study. Although the double-difference method is less sensitive to inaccuracies in the velocity model than standard single event location methods (assuming ray path similarity), there is still a need for consistency in the velocity structure such that for a given event pair, an observed double-difference traveltime residual reflects changes in location rather than velocity. At Decatur, although the microseismicity does not follow a gradual migration away from either injection interval, the timing of events within each region (A-G) is still relatively consistent. This means that the velocity structure will also be relatively consistent during the microseismicity within a given region, ensuring the double-difference residuals can be correctly attributed to their spatial offset. Furthermore, it is also worth considering the extent of the CO 2 plume when assessing its effect on the velocity structure and how this may affect the double difference observations. Strandli et al. (2014) used multilevel pressure data to simulate the CO 2 plume migration, showing that the vertical extent of the plume was limited in depth to approximately 30 m. For two subvertical ray paths from a pair of close events to a common station, there is unlikely to be any significant changes in the traveltime differences if only one of these rays happens to travel through the CO 2 . However, for subhorizontal ray paths, as observed in Fig. 23, then changes to the velocity structure not accounted for in the velocity model may adversely influence the relative locations derived from the double-difference method.
Evaluating the relocation results as a whole, we observe the most significant changes in the lateral position of the event locations. However, significant changes in depth are observed as well, particularly with respect to the removal of the artificial layering of seismicity that was visible in region A prior to the relocation. Moreover, in region C, there is a separation in depth of two distinct subclusters, with the larger showing a sub-vertical lineation previously indistinguishable. In addition, the vertical cross-sections in Figs 9 and 10 show that the seismicity appears much more confined to the base of the Lower Mount Simon and layers below. A possible cause is the lack of faulting higher up in the Mount Simon due to the basement faults terminating within the softer sandstone layers above.
The lateral changes to locations are perhaps most clear in Fig. 11, where the relocation vectors are shown together with the final event locations. Regions A, D, E, F and G all constitute examples of improved alignment in a NE-SW direction, suggesting structural features consistent with this orientation. In making this assessment, it is important to compare these features with the derived location uncertainties in Fig. 20 to ensure any perceived orientation of the structures is not overly influenced by the uncertainties. In the case of region A, a large change in orientation is obtained despite the uncertainties, which favour the prior, more northward orientation. Thus, the true orientation may even be still closer to a more eastward direction. With the exception of a few subclusters in region D, where the uncertainties are smallest, the largest component of uncertainty is not aligned with the observed NE-SW structures in any of the listed regions.
With the seismicity of each region generally becoming less diffuse, with improved clustering and in some cases becoming more linear in orientation, these results will enable an improved interpretation of the pressure interaction with the pre-existing fractures in the subsurface. Using the original, more diffuse locations, Goertz-Allmann et al. (2017a) were able to observe distinct migration patterns of the seismicity in clearly separated regional clusters, suggesting interaction between the lowermost reservoir and basement, possibly caused by lateral heterogeneity in the permeability of the fractured basement. With the regional clusters becoming increasingly discrete, with no clear structural features linking the regions together, our results appear consistent with this interpretation.
Within individual clusters, we now observe much clearer structural features. In region G for example, the 'S'-shaped structure indicates a fault that is not directly observed from the active seismic data. Instead, Williams-Stroud et al. (2020) show two of the faults interpreted from the active seismic data intersecting the region G seismicity, but there is no obvious association between these faults and the occurrence of seismicity. We further observe new lineations within the region G cluster oblique to the main cluster azimuth, potentially suggesting an en-echelon structure, often observed in strike-slip regimes. These observations highlight the complementary nature of microseismicity to conventional active seismic imaging, especially in cases where the latter does not provide sufficient resolution to identify structural features that are crucial for the safe injection of CO 2 . This is particularly true for basement rocks, which are notoriously difficult to image with active seismic studies.
The majority of the microseismicity has been recorded on only the PS3 and GM1 sensors (see supplementary material for the locations based on the operational status of stations). Although these sensors are technically situated in two different wells, they provide almost identical azimuthal coverage. Despite this, the relocations for each region have proved to be stable with the inclusion of the differential back azimuth observations. While the network geometry still plays an important role in the location uncertainties, as shown in Fig. 20, the implementation of the double-difference method has minimized the errors associated with individual event locations.
While we have presented the relocated seismicity with the original catalogue magnitudes in order to directly compare locations, we acknowledge that with changes in location, there will be inevitable changes to the magnitude. However, these are expected to be small due to minimal changes to the source to receiver distances.
Finally, the analysis of the final residuals for each region, in terms of data type, station observations and individual events, greatly benefits how we judge the performance of the inversion and enables us to identify potential data deficiencies and avoid misinterpretation of the results. Despite variations in the weighting of each data type, we observed that overall each data type is being fit reasonably well to the final relocations, with the largest improvements originating from the cross-correlation derived S-wave observations. Although data were carefully selected, analysis of residuals associated with individual stations highlighted that observations on the five GM2 sensors were poorly fitting for regions D and G. Although these represent very few observations compared with those from other stations, and the picking of arrivals from the GM2 sensors was challenging, it is possible that the poor performance is due to a breakdown in the ray path similarity assumption combined with velocity heterogeneities that are not accounted for in the 1-D velocity model. However, given that the GM2 stations are located at relatively shallow depths and the affected regions are relatively close to the sensors, it is unlikely for there to be significant deviations in ray paths for neighbouring events. Despite this, it is possible that further improvements could be made, at least in terms of data fit for the GM2 stations, with the incorporation of a 3-D velocity model. Analysis of residuals for individual events showed that there was consistently good fitting of the cross-correlation derived data for all events, with the main deterioration of residuals associated with the picked P-wave data and to a lesser extent, the back azimuth observations, both of which were weighted relatively low. Removal of events associated with any deterioration in residuals regardless of the data type or weighting showed no significant change to the distribution of seismicity for any of the regions.

C O N C L U S I O N S
We have developed a modified double-difference relocation algorithm, applicable for both 1-D and 3-D velocity models, that incorporates differential back azimuth observations to help constrain the lateral locations of events. This modification is critical for applying the double-difference method to seismicity recorded on a network with poor azimuthal coverage, which is common for the monitoring of microseismicity at fluid injection sites.
Using a 1-D velocity model, we applied our methodology to the microseismicity recorded at the Decatur CCS site, where we relocated 4293 events from an original catalogue of 5397 events. The event relocation involved the manual picking of 114 530 Pwave arrivals, 124 446 S-wave arrivals and 13 902 back azimuth observations. Together with the cross-correlation derived data, the relocation included over 59 million observations for 757 285 event pairs, split across seven geographic regions. Despite the majority of observations being made from the CCS1 and GM1 wells, with an almost identical azimuthal coverage, the results have shown to be reliable with significant residual reductions observed and low uncertainties associated with the final locations. For each region, the relocated seismicity has become less diffuse with improved clustering, and with newly visible linear features often orientated in a near NE-SW direction. We expect these relocations to provide a sound basis for further analysis of the microseismicity triggered by the CO 2 injection, the fluid-subsurface interaction and the interpretation of fractures within the region.

A C K N O W L E D G E M E N T S
We are grateful for support through the Norwegian Climit programme of Gassnova Project Number 618233 and the US-Norway bilateral cooperation agreement on CCS (www.us-norway-ccus.c om). Data are provided by the MGSC. MGSC is funded by the U.S. Department of Energy National Energy Technology Laboratory (NETL) via the Regional Carbon Sequestration Partnerships Initiative (contract number DE-FC26-05NT42588) and by a cost share agreement with the Illinois Department of Commerce and Economic Opportunity, Office of Coal Development through the Illinois Clean Coal Institute. This text has been enhanced by careful reviews from Bob Bauer, Sherilyn Williams-Stroud and Sallie Greenberg of Illinois State Geological Survey, who have proved valuable collaboration partners and provided vital information and feedback on the use of the Decatur microseismic data, the original event catalogue and interpretation of the subsurface.

DATA AVA I L A B I L I T Y
The data underlying this paper were provided by the MGSC as part of a US-Norway bilateral cooperation agreement on CCS. Data will be available from late 2021 via the National Energy Technology

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. Summary statistics for the (a) P-wave and (b) S-wave cross-correlation coefficients for the 4320 events for which reliable observations could be obtained, split into each of the seven regions (A-G). For each box plot, the boxes represent the interquartile range (IQR) for a given region, the blue line within each box represents the median value, while the whiskers represent the upper and lower data limits within 1.5×IQR of the upper and lower quartiles, respectively. Data beyond the whiskers are considered outliers and are plotted as individual blue dots. Figure S2. Summary statistics for the (a) P-wave and (b) S-wave cross-correlation time-lags for the 4320 events for which reliable observations could be obtained, split into each of the seven regions. For each box plot, the boxes represent the interquartile range (IQR) for a given region, the blue line within each box represents the median value, while the whiskers represent the upper and lower data limits within 1.5×IQR of the upper and lower quartiles, respectively. Data beyond the whiskers are considered outliers and are plotted as individual blue dots. A maximum time-lag of ±0.02 s was imposed on the cross-correlation results. Figure S3. Number of observations used in the double-difference relocation for each event pair, split into the five data types for Region A. Figure S4. Number of observations used in the double-difference relocation for each event pair, split into the five data types for Region B. Figure S5. Number of observations used in the double-difference relocation for each event pair, split into the five data types for Region C. Figure S6. Number of observations used in the double-difference relocation for each event pair, split into the five data types for Region D. Figure S7. Number of observations used in the double-difference relocation for each event pair, split into the five data types for Region E. Figure S8. Number of observations used in the double-difference relocation for each event pair, split into the five data types for Region F. Figure S9. Number of observations used in the double-difference relocation for each event pair, split into the five data types for Region G. Figure S10. Map view showing the final set of event relocations coloured by region. The 326 events associated with residuals that show any deterioration from any one of the five data types have been marked in grey. Figure S11. North-south cross-section showing the final set of event relocations coloured by region. The 326 events associated with residuals that show any deterioration from any one of the five data types have been marked in grey. Figure S12. East-west cross-section showing the final set of event relocations coloured by region. The 326 events associated with residuals that show any deterioration from any one of the five data types have been marked in grey. Figure S13. Map views of the final set of event relocations coloured by region and based on the operational status of the stations. Each subplot shows only the events that occurred when the labelled stations were operational. Most events occurred when only the GM1 and PS3 sensors were active, that is Fig. (b).
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.

A P P E N D I X : A N OV E RV I E W O F T H E O R I G I N A L D O U B L E -D I F F E R E N C E M E T H O D
In a linearized inversion for the location of individual events (e.g. Geiger 1910Geiger , 1912, the traveltime residuals r, for an event i,