Full-waveform event location and moment tensor inversion for induced seismicity

Many researchers and industry practitioners have clearly demonstrated that seismic imaging below complex overburden such as salt, for example, requires several key elements that include sufficient target illumination, an accurate velocity model, and an appropriate imaging algorithm (Regone, 2007; Etgen et al., 2009; Gerritsen et al., 2011; Jones and Davison, 2014). Traditional ray-based imaging methods have been superseded by wavefield imaging algorithms, which can inherently handle the multivalued arrivals and wavefront healing expected below rugose and high impedance overburden (Etgen et al., 2009; Jones, 2014). This step change in capability has been augmented by the acquisition of wide azimuth data, which provides improved target illumination and more accurate velocity models. In the case of earthquake event location and source characterization, these same key elements are also required, and one could argue that their inclusion is even more essential if we are to minimize the impact of the earthquakes on society. To address this problem, we have developed two elastic full-waveform workflows, which when used together, provide a powerful solution for subsalt earthquake identification. These workflows are demonstrated by application to shallow and deep monitoring networks, within the giant Groningen gas field in the northeast of the Netherlands.

The Groningen field is the biggest gas field in Western Europe (van Thienen-Visser and Breunese, 2015). The field is operated by Nederlandse Aardolie Maatschappij B.V. (NAM) as a swing producer with production changes made to meet demand. Since 1986, seismic events have been recorded over the north of the Netherlands in what was previously thought to be an inactive tectonic region (Figure ¹). The first seismic event detected near the producing Groningen wells was recorded in 1991 (NAM, 2015). Several medium-magnitude events have been recorded in the region in subsequent years including the Roswinkel earthquake (⁠$MW=3.4$⁠) in 1997 and the Westeremden earthquake in 2006 (⁠$MW=3.4$⁠). However, it was the Huizinge event in 2012 (⁠$MW=3.6$⁠) that received widespread attention. This event was more intense and had a longer duration than previously recorded earthquakes, which resulted in significantly more building damage across the area. The structural damage was aided by the shallow source depth at approximately 3 km and soft surface soils. Multidisciplinary studies, initiated by the Dutch Ministry of Economic Affairs, sought to better understand the relationship between earthquake activity and gas production. These studies concluded that the earthquake activity in the area was a direct result of compaction-induced stress changes due to pressure depletion from production (Bourne et al., 2014).

Since the Huizinge event, seismic monitoring was increased across the field with the installation of an expanded regional, shallow-borehole network in collaboration with the Royal Dutch Meteorological Institute (KNMI), who are responsible for monitoring seismicity in The Netherlands. Although the monitoring network showed a gradual increase in seismic activity, the magnitudes of the earthquakes recorded were relatively small (Figure ¹).

The structural complexity of the subsurface in and around the Rotliegend reservoir plays a significant role in the timing and amplitude of various seismic wave modes propagating from the hypocenter to the receivers. Figure ² displays a cross section through the subsurface velocity model, which highlights the structural complexity. The subsurface geology is broadly divided into the following main units: the suprasalt section, which includes the North Sea formation of clastics overlain onto Cretaceous carbonates, followed by clastic Jurassic and Triassic formations; the Zechstein evaporite section, which includes a variable-thickness floating anhydrite layer and basal anhydrite (average thickness of approximately 50 m); and the Rotliegend reservoir interval comprising mainly Aeolian sandstones toward the top with more fluvial deposits of sand, conglomerates, and clay toward the bottom. The porosity within the reservoir can be up to 25%, and it is a contributing factor for induced seismicity because areas of high porosity are likely to compact and subside following gas extraction (NAM, 2013), which may increase vertical stress levels (Ketelaar, 2009) making the chance of inducing an earthquake more likely.

Below the reservoir lies the Carboniferous, composed of mainly red bed facies of sandstones and shales (Stauble and Milius, 1970). From a field development perspective, the basal anhydrite and Zechstein salt are the seal to the Rotliegend reservoir, which has been charged from the coal-bearing Carboniferous below (Stauble and Milius, 1970).

A key requirement in the study of the induced-seismicity problem is to accurately locate earthquakes. There are many methods that have been developed for earthquake location. These include the hypocenter method (Aki and Richards, 1980); the double-difference method (Waldhauser and Ellsworth, 2000; Zhang and Thurber, 2006); the differential time method, also referred to as the equal differential time (EDT) method (Font et al., 2004; Lomax, 2005; Satriano et al., 2008; Theunissen et al., 2012; Spetzler and Dost, 2017); and the grid search method (Rodi, 2006; Pickering, 2014). All these methods either use traveltime differences between stations to estimate location, or search over a prescribed grid of potential earthquake locations and use minimization criteria to determine the most likely location.

A collapsing grid search method was used by Pickering (2014) to derive locations at the Groningen field using events from two deep boreholes that penetrate the reservoir. In this method, raypaths were calculated as a vector in terms of distance, azimuth, and dip. Instead of iterating linearly toward a minimum residual solution, the grid search algorithm searches over a coarse grid for a point that minimizes the misfit between the picked traveltimes and the theoretical traveltimes produced by ray tracing through the subsurface velocity model. This is applied for every receiver and for every source location. When a source position is found that minimizes the sum of the square residuals, a smaller, more finely sampled grid volume is defined around the initial source location. It is assumed for each search grid, that the located source point in the subvolume is the closest grid point to the global minimum. The algorithm continues in this fashion until some user-specified location resolution is obtained. Although the application of this method can potentially offer more accurate depth resolution, depending on the spatial sampling of the modeled grid, it has the same limitations as previously noted, as applied by Pickering (2014), i.e., the use of a 1D laterally invariant velocity model and the choice of single traveltime arrivals. Grid search methods are also typically computationally intensive compared with iterative solutions based on equation 2.

More recently, Spetzler and Dost (2017) apply the EDT method to determine locations of events recorded by the shallow borehole network. In the EDT method, the origin time of the earthquake can be determined after the focus of the earthquake has been estimated. The recorded arrival time for a wave to travel from the source to the receiver is estimated by integrating along a raypath to produce the recorded arrival time. The unknown origin time is eliminated by subtracting the traveltime from two receiver positions for the same source location. In their paper, Spetzler and Dost (2017) solve a depth-dependent objective function, which uses local 1D velocity profiles to compute a series of traveltime functions, for a range of source depths and epicentral distances. A single-arrival eikonal solver is used to calculate the traveltimes. An advantage of the authors’ approach is that a total misfit function is obtained by summing contributions of differential traveltime residuals for all station pairs. However, the reliance on 1D invariant models and single-arrival traveltime (direct P-wave only) makes the method insufficient for the Groningen subsurface problem, as will be demonstrated later.

To overcome the restriction of laterally invariant velocities and methods that use single-event arrivals or that require a minimum of three active boreholes, we have adopted full-waveform methods and applied these to the Groningen borehole networks. The term “full waveform” can mean different things, for example, the use of finite-difference schemes, using different components of the wavefield, or using elastic properties of the subsurface. Migration-based methods, for example, have sought to focus the earthquake energy by back propagation of the recorded wavefield in time or space to obtain a focused event location (McMechan, 1982; Rentsch et al., 2007; Lu and Willis, 2008; Xuan and Sava, 2010). These methods usually require a densely sampled geophone network to sufficiently obtain a focus and minimize migration artifacts (Bazargani and Snieder, 2013). Alternatively, waveform inversion-based methods (Song and Toksöz, 2011; Li, 2013; Li et al., 2016) can be used to jointly solve for the location and source mechanism by waveform matching and moment tensor inversion.

In this paper, we use the methodology developed by Li (2013), but we apply this in separate workflows that are based on either waveform attributes (arrival time and P polarization) or the waveforms themselves. We use the elastic properties of the subsurface, multicomponent recordings and apply these in a finite-difference waveform inversion scheme. We show that, in fact, the subsurface wave propagation at Groningen is very complex and that this requires 3D elastic traveltime estimation, which uses all arrivals to improve location accuracy. Furthermore, the seismic focal mechanism is an integral element of earthquake characterization, and it can also be derived simultaneously along with location.

To sufficiently detect seismic activity, several monitoring networks have been installed across the Groningen field (Figure ³). These networks provide the opportunity to detect earthquakes from within the reservoir via two deep boreholes and over the full lateral extent of the field via a shallow borehole network.

The deep borehole network is composed of two deep (3000 m) boreholes that are instrumented with multicomponent (3C), 15 Hz geophones that penetrate the reservoir to approximately 3 km depth (Figure ³). These wells were converted from old monitoring wells near the towns of Stedum (SDM-1) and Zeerijp (ZRP-1). The SDM-1 well was equipped with ten 3C geophones, and ZRP-1 was equipped with a string of seven 3C geophones. Sensors in the string are positioned 30 m apart, spanning a total length of 270 and 180 m, for SDM-1 and ZRP-1, respectively. In SDM-1, the geophone depths are at 2755–3025 m (below mean sea level [MSL]), with the first geophone located within the basal anhydrite layer and the other geophones positioned in the Rotliegend reservoir section (see Figure ³). In ZRP-1, the geophone depths are at 2785–2965 m (below MSL), with the first two geophones within the basal anhydrite layer, and the other geophones are in the reservoir section (Figure ³). For all practical purposes, the ZRP-1 well can be regarded as a vertical well, the SDM-1 well is slightly deviated in the monitoring section, which must be considered when the orientations of the geophones are derived. The geophone strings use Sercel Slimwave tools, equipped with OMNI-2400 3C, 15 Hz geophones measuring orthogonal $x$⁠, $y$⁠, and $z$ components in a left-handed coordinate system.

At the surface, a shallow (200 m depth) regional borehole network has been installed by NAM comprised 70 boreholes on a 6 km staggered grid, which provides an interborehole spacing of 3–4 km (Figure ³). The boreholes are instrumented with four 3C 4.5 Hz geophones down the wells and a single accelerometer placed at the surface. The network has been handed over to KNMI, who are responsible for monitoring the seismic activity. The geophones in the string are spaced 50 m apart as shown in Figure ³.

These networks represent one of the most detailed seismic monitoring efforts ever undertaken over a producing field. The recorded events have been studied to determine the microearthquake event locations and focal mechanisms.

The velocity model used for this study has been developed over many years using prestack tomographic inversion, and it has been constrained with sonic velocities from the multitude of boreholes that have been drilled in the region (Romijn, 2007). The seismic P-wave velocities within the reservoir section are higher at the top and bottom and decrease toward the middle (⁠$3800m/s$ average); in comparison, the overlying anhydrite has a high velocity of $5900m/s$ (Figure ² and ²). In addition, the average upper Carboniferous velocity immediately below the reservoir is $4250m/s$⁠. The lower velocity reservoir is, therefore, bounded by high impedance contrasts at the top and bottom, which effectively channels earthquake energy within the reservoir section and causes significant mode conversion as the seismic energy radiates away from the hypocenter. At shallower depths, velocity contrasts are present between the overburden and the Zechstein salt, as well as contrasts within the suprasalt section, including across the Cretaceous chalk interval and between the upper and lower North Sea formations. This all results in complex wavefield propagation, which needs to be correctly modeled and interpreted before a complete understanding of the earthquake dynamics can be achieved.

To investigate the wave-propagation phenomena in more detail, Figure ⁴–⁴ displays the wavefront snapshots from 3D elastic, finite-difference simulations using a double-couple source and for different positions across the reservoir interval. The modeling software is a proprietary package and uses eighth-order accuracy in space and second-order accuracy in time. The simulation was run up to a maximum frequency of 20 Hz using a Gaussian wavelet. The examples shown simulate the propagation for an induced earthquake from the reservoir to the surface, as it would be recorded by the shallow borehole network. The recorded wavefronts contain P- and S-arrivals as well as significant mode-converted arrivals with P-S and S-P conversions. This is clearly seen in Figure ⁴–⁴, which includes event gathers recorded along the surface for a duration of 4 s from the earthquake initiation. The mode-converted events can be traced to the top Zechstein salt interface and across the Cretaceous chalk boundary with the overlying North Sea formation. Several different arrivals are therefore present within the seismic record. If multiples are also included, then the seismic recordings are very difficult to interpret, as KNMI have found (Spetzler and Dost, 2017). The simulations show that a significant proportion of the earthquake energy travels horizontally along the reservoir and that mode conversion is largely controlled by the impedance contrast at the top and base Zechstein, anhydrite floaters, as well as chalk boundaries. As previously mentioned, the variable thickness of the top salt as shown in Figure ³ is also an influencing factor and causes strong wavefront divergence.

In the case of the downhole monitoring stations, one can expect to record several different refracted and reflected arrivals within the reservoir section, as displayed in the elastic wavefield simulation of Figure ⁵ and the seismic trace (Figure ⁵). For this simulation, we increase the maximum frequency to 100 Hz to show the detail in the wavelet. The first arrival in this case is a weak P head wave transmitted along the reservoir and overlying anhydrite. This is followed by a weak P-wave direct arrival and then by a very strong guided P-wave, a direct S-wave, and finally a very strong guided S-wave. The presence of guided waves makes direct P- and S-wave arrival interpretation very challenging. The significantly higher amplitude of the guided waves can be erroneously interpreted as direct arrivals, which would potentially result in a location error. As previously mentioned, single arrival event location algorithms would struggle in this situation and event locations, for example, derived from the EDT method, would be unreliable with only two wells available to constrain the solution.

An underlying observation that should be considered when locating events is that there is a high dependence on where one observes the earthquake, considering that the network may not fully capture the radiation pattern of the event. The choice of surface or downhole monitoring is a key consideration, but also, because the subsurface is strongly heterogeneous, the source-to-receiver azimuth and offset are important factors. This will be touched on later in the discussion of error estimation.

To overcome the complex wave phenomena highlighted by the elastic modeling simulations, elastic-wave-based event location workflows, which use all arrivals of the recorded data, have been applied for the Groningen events.

It would be desirable if waveforms can be directly used to locate events and, preferably, to derive the source mechanism via moment tensor inversion. However, this is not an easy task (Li, 2013). First, if observed waveforms are used directly, we need to generate synthetic waveforms that have all the complexity as observed in the field before a direct comparison can be made. Because we do not know the source locations upfront, we need to generate synthetic traces from all possible source locations (i.e., grid search). Second, because our velocity model is only an approximation to the true subsurface structure, the higher the frequencies we want to use and the larger the source-receiver distance, the greater the phase mismatch between the observed and synthetic data could become. Therefore, without appropriately conditioning the data, any direct waveform matching algorithm is unlikely to succeed. To overcome this problem, we apply two methodologies in a cascaded fashion one after the other and we avoid using high frequencies in the inversions. The first method (the waveform attributes method [WAM]) uses the kinematic information retrieved from the synthetic and observed waveforms to provide a preliminary event location, whereas the second method (the full-waveform method [FWM]) automatically determines the source parameters by matching the synthetic waveforms with the observed ones to generate the final location and to provide an event focal mechanism via moment tensor inversion.

The basic concept of the wave-equation event location method is based on a grid search in space (and potentially additional source attributes), which minimizes the attribute difference between the observed and synthetic data. In the WAM approach, the minimization is undertaken on waveform attributes derived and extracted from the data, more specifically, the P- and S-wave traveltime and the P-wave polarization. A full description and mathematical foundation of the methods can be found in Li (2013), and it is summarized in the following paragraphs.

The overall workflow is schematically shown in Figure ⁶. The input field data are prepared by applying tool rotation to orient the receivers to the incident wavefield and converting to the northeast down geographical coordinate system. Polarization analysis of check-shot data was applied to derive the instrument rotation angles. Three check shots were used to orient the instruments for the SDM-1 and ZRP-1 wells. Some preprocessing of the field data is applied to remove unwanted noise by band-pass filtering (using passbands of 1–10 Hz for the shallow network and 8–25 Hz for the deep network), random noise attenuation and to compensate for instrument roll-off by applying instrument amplitude and phase corrections. The P- and S-arrivals have been automatically picked but may be adjusted manually in areas with poor S/N. Next, synthetic data are generated by using 3D elastic finite-difference forward modeling for the Green’s function calculations over a dense grid (⁠$50×50×25m$⁠) of candidate source locations. The use of full 3D elastic modeling is applied, so that all wave modes can be used in the waveform solution. To reduce the run time for the elastic wave simulations, the calculations are limited to a $12×12km$ box around each well location for the shallow borehole data set and a $5×5km$ output area for the downhole data set. The vertical extent of the source grid is limited to an area above and below the reservoir (2500–3500 m). Higher frequencies are modeled, compared with the inverted range of the data, to study the synthetic waveforms at large distances. A maximum frequency of 100 Hz is used for the downhole simulations, whereas a lower frequency wavelet up to 20 Hz is used for the surface network modeling. The run time is also optimized by running the modeling in reciprocal mode thus reducing the number of modeled earthquake locations. A comparison of the observed and modeled data is shown in Figure ⁷ for the three geophone components of the bottom receiver (200 m depth) from eight shallow boreholes located around the earthquake epicenter. In this comparison, the P and S picks are also shown as indicated by the red and blue picks. The synthetic traces have been generated with an absorbing upper boundary and therefore do not include surface multiples.

The event location portion of the workflow incorporating the $F1$ and $F3$ terms previously discussed is then applied, and the global minimum is selected, which represents the final event location in the $x$⁠, $y$⁠, and $z$ spatial coordinates. A comparison of the various contributions of the $F1$ and $F3$ terms for the two downhole wells previously discussed are shown in Figure ⁸. The $F1$ term, which represents the difference in traveltime between the observed and synthetic data, forms a concentric minimum around the well locations. Using traveltime differences alone would be insufficient to locate a global minimum unless the network contained at least three wells, as was previously determined by Vavrycuk (2007). The polarization term $F3$ indicates that the P-wave incident wavefield and the objective function have two unique solutions, which are separated by 180°. As with the $F1$ term, additional well contributions would improve the location accuracy when using polarization alone; however, this is also nonunique for the two wells’ geometry. Adding the $F1$ and $F3$ terms for both wells provides a unique solution as shown in Figure ⁸ (bottom right panel). In practice, we use the WAM workflow to determine a preliminary event location, which is then used to reduce the search area for the FWM workflow, which is now described.

The full-waveform event location workflow that we use was introduced by Li (2013) and Li et al. (2016). In this method, a synthesized, fully elastic earthquake waveform is matched in a least-squares sense to the waveform recorded in the field at the geophone locations. The end-to-end workflow is schematically shown in Figure ⁶. As with the WAM workflow, the input data are prepared in the same way.

The waveform inversion uses an alignment and matching algorithm, which removes residual misalignment between synthetic and observed traces using an optimized crosscorrelation method. This algorithm produces reference double-couple solutions as regularization terms for the subsequent inversion (Li et al., 2016). Therefore, an advantage of the FWM workflow is that it only requires approximate P- and S-arrival times. In practice, we apply a single window for the alignment, which encapsulates the P- and S-arrivals; this ensures that the P-to-S separation is preserved.

After the $T0init$ correction has been applied, the remaining shift for alignment $Δτn$ for each trace $n$ can be found by crosscorrelating a very large number of trial solutions with the observed data following the procedure described in the fast alignment and matching algorithm of Li et al. (2016) and is summarized here:

1. 1)

   For each potential event location $(x,y,z)$⁠, crosscorrelate the observed data with the synthetic traces for each receiver. This essentially involves crosscorrelating tens of thousands of candidate solutions with the observed data. We assume a double-couple focal mechanism, and we generate trial solutions for every 5° in strike, dip, and rake dimensions.

2. 2)

   For each trial solution, find the maximum of the normalized crosscorrelation series and the corresponding residual delay. The residual delay is subject to a constraint to prevent errors and is described in the next paragraph.

3. 3)

   Stack the maximum normalized crosscorrelation values for each trial solution.

4. 4)

   The maximum normalized crosscorrelation values are then compared for all locations and receivers, and the global minimum is selected. Correspondingly, the best estimated time delay is also selected.

To expedite the automatic search over all candidate source locations, we evaluate the root-mean-square (rms) error between picks from observed and synthetic data, and if it is larger than a given threshold, we consider this position unlikely to be the correct source location, save the evaluation, and move on to the next candidate position. Note that any misalignment between the two wavelets may arise from two potential sources. The first cause is the rough estimation of the origin time $T0init$ from traveltime picks only. This achieves bulk alignment of the traces for the crosscorrelation, but it also inevitably introduces some misalignment at the same time. The second cause is from inaccuracies that may be present in the subsurface velocity model, as previously mentioned.

It is well known that full moment tensor inversion is often ill conditioned for sparse observation geometries (Rodi, 2006; Vavrycuk, 2007; Jechumtálová and Eisner, 2008; Eaton and Forouhideh, 2011; Song and Toksöz, 2011). We therefore prefer to include the stabilization factor as mentioned in equation 10 to better control the performance in these conditions. An example of the objective functions for these terms is shown in Figure ¹⁰ for the two deep boreholes. As can be seen from the images, the objective function over the depth slice is elongated in a direction that is orthogonal to the azimuth between the two wells. The individual terms for the objective function contribute to locating the event hypocenter, but it is the summation of all terms that provides the final location as shown in the bottom right image of Figure ¹⁰. In this example, the contribution from the $RM$ term is very weak and difficult to interpret because there are only two wells available to sample the earthquake’s radiation pattern. Several different interpretations are, therefore, possible for this focal mechanism, and it is necessary to check the waveform alignment carefully for consistency around the final location.

The final stage of the workflow is to fully quality control the event location and moment tensor inversion results. This is achieved via a custom-built interactive software tool that allows one to investigate the objective function space and assess the data misfit by comparing the observed and synthetic displacement traces for any event in the database. Interactivity is further provided to allow for alternative scenario testing of either the event location or focal mechanism. This provides a powerful diagnostic capability that helps the user to rapidly eliminate any improbable results and to determine the best fit to all the available data. At Groningen field, we are expecting the source mechanisms to indicate dip-slip, normal motion, as the reservoir compacts from production of the gas (NAM, 2015).

As previously mentioned, we use the WAM workflow to generate an initial event location and the FWM method to provide a final location and focal mechanism. If the S/N of the data is good, it is possible that the WAM method may be sufficient. However, as has been demonstrated, the wavefront arrivals that are encountered in the Groningen subsurface are complex, and this requires a more complete use of the waveform coda than just P- and S-arrivals alone. To demonstrate the added benefit, Figure ¹¹ compares the location accuracy for the WAM and FWM workflows against a ray-based solution for a known event. A known event in this context refers to a synthetic event generated from the Green’s functions and a specified moment tensor. As can be seen in Figure ¹¹, the ray-based solution, which was generated with ray tracing rather than elastic wave modeling, shows an event that is mislocated from the true location. Application of the WAM workflow produces a location, which is closer to the correct position, i.e., closer to the fault, but the FWM workflow generates an event that is colocated with the true position. The horizontal location difference between the ray-based and FWM event position is approximately 350 m and varies by 175 m in depth.

The deep borehole network includes the SDM-1 and ZRP-1 wells (Figure ³). The event catalog for these wells spans the recording period from 10 October 2013 to 29 June 2015, and it contains a total of 522 events. Due to instrument failures, the geophone strings were pulled from the boreholes, repaired, and redeployed, which resulted in different recording phases for SDM-1 and ZRP-1. Although it is not a strict requirement for the event location portion of the workflow, we selected data that were recorded by both wells simultaneously to improve the location accuracy and constrain the moment tensor inversion. In addition, events with good S/N and magnitudes between $−0.7$ and 1.5 were selected for the location workflow. This brought the final number down to 45 events that were located (Appendix  A). The results are shown in Figure ¹² along with the focal mechanisms. The results show good positional alignment to the major fault trends that offset the reservoir. However, it is clear from the derived focal mechanisms that there is a broad range of dip and strike orientations that do not appear to be consistent or representative of what we would suspect from normal faulting. This is verified by plotting all the results together on a stereonet display (Figure ¹³). Here, we used the freely available software developed by Cardozo and Allmendinger (2013) for the stereonet displays and plot the density of all the derived dip/strike values from the moment tensor inversion (Figure ¹³) and for the primary and auxiliary planes (Figure ¹³). There are no obvious orientations that match the fault trends shown in Figure ¹². To explain these results, we more closely investigated the objective functions used to constrain the event location and the focal mechanism. An event was selected in the center of the survey for this purpose (Figure ¹⁴). The event location (Figure ¹⁴) is at the minimum of the objective function; however, the shape of the objective is elongated orthogonal to the azimuth between the two wells. Essentially, this area covers the overlap between each well and effectively defines the axis along which we would expect the largest location errors to occur.

Analysis of the focal angles, as used in the source mechanism objective function (Figure ¹⁴), indicates a multimodal objective function. The selected minimum in this example has fallen in a region of the focal angle domain (strike = 216, dip = 26, and rake = 120°) that indicates thrust faulting, which is not expected for the Groningen field. This result was generated with no constraints; i.e., there were no restrictions placed on focal angles for the dip, strike, or rake angles. If normal faulting is assumed, then it would be relatively easy to constrain the solution to only select minima over a prescribed angle range. Although this is an option, if not done carefully, it might bias the result to what we want to see rather than what the data are suggesting. Of some concern is the multiminima seen in the objective function. We hypothesize that this is due to the use of only two wells that have limited source-to-receiver azimuth coverage and cannot fully capture the radiation pattern of the earthquake. Essentially, the earthquake focal sphere is not sampled well enough, so the true characteristics of the earthquake cannot be determined, and the solution is nonunique. To resolve this ambiguity, we must monitor the earthquake at more locations as suggested by Rodi (2006), Vavrycuk (2007), Jechumtálová and Eisner (2008), Bardainne et al. (2016), and Willacy et al. (2018).

The results from 100 located events using the full-waveform workflow (Appendix  B) are displayed in Figure ¹⁵. These events were recorded by the shallow borehole network between 2015 and 2017, and they are shown superimposed on the top reservoir (Rotliegend) depth map, which shows the location of the major faults across the area. The events were selected based on good S/N and cover a range of magnitudes, $ML$⁠, from 0.1 to 2.2. For the Green’s function modeling, only the deepest receiver was used (to avoid near-surface noise), with a spatial aperture of 12km and a maximum wavelet frequency of 20 Hz. Using much bigger apertures becomes very costly for 3D elastic simulations, but a 12 km aperture allows up to eight wells to be included in the objective function, which provides a good coverage over a broad range of azimuths.

The results show a strong correlation between the earthquake location and the proximity to the faults (Figure ¹⁵). Most of the events lie on the major north–northwest to south–southeast-trending faults with a few located on smaller faults that appear to join these larger faults at right angles. The source mechanisms, as determined by the moment tensor inversion, are also remarkably consistent and match with the dominant fault trends and show a predominantly normal sense of motion as shown on the stereonet plot in Figure ¹⁶. A few of the events appear to indicate a more strike-slip motion. In some cases, these events are located at the junction with several other faults, and it is uncertain whether the motion takes place on more than one fault simultaneously. This needs further geomechanical evaluation.

As previously undertaken for the deep borehole network, we now take a closer look at one event detected by the shallow borehole network and analyze the objective function for location and source mechanism (Figure ¹⁷). The objective function, which is the summation of contributions from eight shallow boreholes, shows a very well-defined minimum (Figure ¹⁷). Automatic selection of the minimum is easily achieved and provides increased positional accuracy compared with the previous two-well, deep-borehole results.

The objective function for the focal mechanism (Figure ¹⁷) also contains a well-defined minimum. The selected focal mechanism for this example has focal angles of dip = 55°, strike = 170°, and rake = −75°. This is consistent with normal faulting, and the orientation closely matches the strike and dip of the nearest interpreted fault.

The results from the full-waveform workflows provide a consistent interpretation for the events analyzed, where the events have good receiver coverage and corroborate previous studies in the area (Bourne et al., 2014). However, for any methodology, it is important to get a sense of error in the position and/or depth. This can be somewhat subjective because the criterion and method for different algorithms varies as does the perspective of the interpreter to what constitutes an error. We have analyzed various types of errors within the FWM workflow previously described, including errors in velocity, grid sampling, and event traveltimes. Each location will have its associated errors depending on the signal quality, velocity model accuracy, and recording coverage. To judge the overall sensitivity in the spatial position, we have derived an error ellipsoid from the variance of the objective function. This has been achieved by taking a selection of data points from around the minimum location of the full-waveform objective function and measuring the spatial variance and standard deviation. We use a range of 20% of the waveform correlation value to select the data range. Once these points have been selected, the errors are computed by deriving the eigenvalues and eigenvectors of the spatial covariance matrix (Montalbetti and Kanasewich, 1970). Intuitively, the covariance matrix generalizes the notion of variance to multiple dimensions. This provides a quantitative estimate of the variance of the objective function around the minimum. Figure ¹⁸ shows an example depth slice through the full-waveform objective function computed from four shallow wells. The clipping of the objective function based on the correlation value is shown in Figure ¹⁸. As can be seen from the displays, the four wells provide a limited azimuthal range around the event epicenter. The minimum of the objective function forms an ellipsoid from which we can measure the axial variation; in this example, the largest variation is in the northeast to southwest direction. The axis of the error ellipsoid corresponds to the eigenvectors of the covariance matrix and their lengths to the square roots of the eigenvectors. Figure ¹⁸ shows another event location, but in this example, there is good azimuthal coverage from the 10 surrounding wells. The minimum of the objective function is therefore spatially localized, and the variation in location error is small (Figure ¹⁸). The error analysis has been applied to events recorded on the shallow borehole receivers that have good azimuthal source-to-receiver coverage (Appendix  B). The results indicate that the spatial errors are on average <$100m$ with a standard deviation of 50 m.

The dependence on the borehole coverage was also noted by Willacy et al. (2018), who present the location accuracy and focal mechanism sensitivity by analyzing the results of different combinations of input well contributions to the objective functions. The authors hypothesize that the azimuthal coverage is effectively being limited when the number of wells are reduced; therefore, a smaller portion of the earthquake radiation pattern is being recorded. This was demonstrated mathematically by Vavrycuk (2007), who concludes that a minimum of three wells were needed to derive an accurate moment tensor from P-wave amplitudes alone. If P- and S-waves were used, then this could be reduced to two wells. However, we have found that in the case of the Groningen subsurface, the complex wave propagation results in a nonunique solution, so either a priori constraints must be used, as with some approaches (Song and Toksöz, 2011), or additional deep monitoring wells are required.

The proximity to the faults has also been analyzed by reference to the top Rotliegend depth map. We have noted that the events derived from the FWM are statistically on average 55 m closer to the faults than the published KNMI catalog values from the 150 events analyzed. However, it should be noted that this is somewhat subjective and depends on the ability to interpret the fault offsets on the top reservoir surface.

The estimation of the depth error from the objective functions results in an average value of 180 m with a standard deviation of 125 m. Perhaps a better way to visualize the scatter in depth position is to display the depths relative to the reservoir interval. This is undertaken in Figure ¹⁹ for the 150 events detailed in this paper. So far, all events to date have been located within the reservoir section; however, based on the error analysis, we cannot with certainty rule out the possibility that some events could occur outside the reservoir interval. Spetzler and Dost (2017) locate a few events above the reservoir using their location method and postulate that these events may have occurred by brittle rupture of the top anhydrite layer. Although mechanically this is plausible, the depth errors associated with their technique may be too large to provide any degree of certainty. Combining the shallow and deep networks within the same inversion may provide additional positional accuracy. However, this is left for future work.

We have applied a full-waveform event location workflow, which implicitly uses all arrivals and provides for event location and source mechanism simultaneously. The results from analysis of 150 microearthquakes recorded by the deep and shallow borehole networks demonstrate a very close correlation to the major faults, which offset the primary Rotliegend reservoir interval with a north-northwest–south-southeast strike direction. The sense of motion, as determined by moment tensor inversion, is consistent with normal faulting. In comparison, we also see a close correlation to the event locations derived from the deep borehole network installed in the SDM-1 and ZRP-1 wells; however, the limited azimuth coverage of the deep boreholes means that the focal sphere is undersampled, which causes multiple minima in the objective function and the derivation of source mechanisms is nonunique.

By using an elastic wave equation to solve for the event location and source mechanism simultaneously, we have overcome the inherent problem of requiring accurately interpreted arrivals and avoid the single arrival limitation of previously published methods by using all arrivals to improve the accuracy of the located event. This approach is more reliable and able to use poorer quality data because the match is solved in a least-squares sense. However, the methods presented in our study require an accurate velocity model. This has been achieved for the Groningen subsurface due to the many wells that have been drilled since production started and which have been incorporated into the tomographic velocity inversions, as part of the prestack depth imaging workflows for the active source, surface data sets.

All our derived events to date have been located within the reservoir bounds. Although there are bigger vertical errors due to the detection coverage, we think that it is not likely that events will be located outside of the reservoir, although at present this cannot be ruled out. It is hoped that by using the surface and downhole networks simultaneously, we will see an improvement in the vertical positional accuracy. However, from a practical perspective, the shallow and deep borehole networks have limited operational time overlap, so there are relatively few events in which a comparison can be made.

The results from this study support the previously published work, which suggests that reservoir compaction-generated fault motion, via gas production, is the primary cause of the earthquakes. Now that we have established a set of accurate locations, these can be used to constrain the geomechanical models of earthquake rupture. Additional work is ongoing on this topic, and it is hoped that further results will provide more insight into the effects of induced seismicity at Groningen.

The authors would like to thank J. van Elk, D. Doornhof, J. Tomic, and R. Romijn from NAM. We would also like to acknowledge our Shell colleagues F. ten Kroode, S. Oates, and W. Mulder for many productive discussions and guidance throughout this work. Our thanks also go to KNMI for providing the geophone and accelerometer data and to NAM and Shell for providing permission to publish the material presented. This paper benefited from suggestions provided by the anonymous reviewers. The event catalog mentioned in the "Results" section is available at http://rdsa.knmi.nl/dataportal/.

Data associated with this research are confidential and cannot be released.

Appendix A contains the event location results derived from the deep borehole network using the workflows presented in this paper. The events are provided with the KNMI catalogue information, which includes event time and magnitude.

Appendix B contains the event location results derived from the shallow borehole network using the workflows presented in this paper. The events are provided with the KNMI catalogue information, which includes event time and magnitude.