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

Locating microearthquake events below complex heterogeneous overburden requires robust location methodologies that can honor multipathing in the seismic wavefield. We have developed two full-waveform event location methods that form a complementary solution for locating earthquakes and simultaneously deriving focal mechanisms via moment tensor inversion. The methods are based on the application of 3D elastic wavefield modeling, which is used to generate waveforms and extract wavefield attributes, for comparison to the observed field data. Events are located and focal mechanisms are derived via a multiparameter inversion, which minimizes the differences between synthetic and observed data. The results have been applied to the induced seismicity observed within the giant Groningen gas field, onshore Netherlands, where recorded earthquakes are triggered by stress changes, induced in the reservoir through pressure depletion. Locating events below the field is compounded by the presence of strong guided waves, which are trapped in the lower velocity reservoir interval. This complex multivalued wavefield is problematic for traditional event location methods, which assume a single traveltime arrival. We overcome this limitation by using all event arrivals in a wave-based solution to improve the accuracy of locating earthquakes and overcome the ambiguity of solving for location and the focal mechanism simultaneously. The event location methods have been applied to shallow and deep monitoring networks, and 150 events have been located with high accuracy. The interpretation of the earthquake activity indicates that the events studied originate from the movement of larger graben bounding faults, which are oriented in a north-northwest– south-southeast direction.


INTRODUCTION
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 raybased 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 1).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 (M W ¼ 3.4) in 1997 and the Westeremden earthquake in 2006 (M W ¼ 3.4).However, it was the Huizinge event in 2012 (M W ¼ 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 1).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 2 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 inducedseismicity 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.
At Groningen, events recorded by the regional shallow monitoring network (Figure 3a) are transmitted to KNMI, who perform an initial event location using a semiautomatic workflow based upon the hypocenter approach (Aki and Richards, 1980;Lienert et al., 1986).In this method, the calculated traveltime t C i at station i can be written as where T is the traveltime as a function of the location of the station ðx i ; y i ; z i Þ and the hypocenter.This equation has four unknowns, namely, the origin location ðx 0 ; y 0 ; z 0 Þ and time ðt 0 Þ.To solve this equation, recordings from at least three stations are needed to determine the hypocenter and origin time.This is normally solved by calculating the misfit or residual r i between the observed and calculated traveltimes, which is equivalent to the difference between the observed and calculated arrival times where t 0 i is the observed arrival time.In principle, this is a simple task; however, the traveltime function is nonlinear and must be calculated.In the case of KNMI's workflow, this is achieved by ray tracing through an average 1D velocity model (Spetzler and Dost, 2017).The success of this method, therefore, depends on the accuracy of the traveltime estimation and the ability to accurately interpret the arrival types, i.e., P-and S-waves. Figure 3a, also shows the Zechstein salt isopach map, which strongly influences the amount of ray bending, due to the velocity contrast between the salt and surrounding rock.For the shallow borehole network, as will be discussed later, the interpretation of multiple arrivals is, therefore, not straightforward, and KNMI elect to use the direct P arrival only as the direct S-wave is difficult or impossible to interpret by hand.In addition, in the hypocenter method, the depth of the event is required to calculate the traveltime, and in practice, this was fixed at an average reservoir depth of 3 km for the Groningen area (Spetzler and Dost, 2017).Knowing the exact hypocenter as well as the epicentral location is required to enable an accurate geomechanical interpretation of the field.
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 singlearrival traveltime (direct P-wave only) makes the method insufficient for the Groningen subsurface problem, as will be demonstrated later.

Full-waveform event location
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.

MONITORING INDUCED SEISMICITY
To sufficiently detect seismic activity, several monitoring networks have been installed across the Groningen field (Figure 3a).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 3b).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 3b).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 3b).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 3a).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 3c.
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.

ELASTIC WAVEFIELD MODELING
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 (3800 m∕s average); in comparison, the overlying anhydrite has a high velocity of 5900 m∕s (Figure 2a and 2b).In addition, the average upper Carboniferous velocity immediately below the reservoir is 4250 m∕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 4a-4c 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 Pand S-arrivals as well as significant mode-converted arrivals with P-S and S-P conversions.This is clearly seen in Figure 4d-4f, 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 3a 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 5a and the seismic trace (Figure 5b).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.

FULL WAVE EVENT LOCATION METHODOLOGY
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.

WAM
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 Pwave 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 WAM objective function, Fðx; y; zÞ, to be minimized is defined using two terms of Li's original formulation, namely, F 1 describing the minimization of the observed and synthetic traveltimes and F 3 describing the minimization of the polarization for the observed and synthetic data.These terms are calculated as follows: Fðx; y; zÞ ¼ λF 1 þ ð1 − λÞF 3 ; (3) where (4) in which the difference in moveout between observed and synthetic (P and S) traveltimes is minimized and T init 0 is calculated to align the observed and synthetic arrival times (this can be calculated per well or for all wells simultaneously; h : : : i; depicts the average): and the second term, which calculates the square root of the summation over all stations i of one minus the inner product between the P-wave polarization unit vectors of the observed data (e obs P;i ) and synthetic data (e syn P;i ) at station i.The parameter λ in equation 1 specifies the weight assigned to the traveltime difference term F 1 , and it can be chosen between zero and one.The choice of weight should consider the Figure 6.Full-waveform event location workflows used in our study.(a) The WAM and (b) the FWM.The steps highlighted in gray are different between the two workflows.In practice, we use the WAM to determine an initial event location and then use this position as a starting location to refine the search box for the FWM.quality of the polarization, which may be ambiguous at large offsets.For the data studied, we used a value of λ ¼ 0.5.The misfit function is evaluated at every grid point of the volume defined by the Green's function library (synthetic hypocenter data grid).
The overall workflow is schematically shown in Figure 6a.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 × 25 m) 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 × 12 km 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 7 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 F 1 and F 3 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 F 1 and F 3 terms for the two downhole wells previously discussed are shown in Figure 8.The F 1 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 F 3 indicates that the P-wave incident wavefield and the objective function have two unique solutions, which are separated by 180°.As with the F 1 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 F 1 and F 3 terms for both wells provides a unique solution as shown in Figure 8 (bottom right panel).In practice, we use Figure 7.A comparison of observed and synthetic displacement data for the three geophone components extracted from eight wells around the epicenter location and sorted by offset distance.The P-and S-arrivals are shown by the red and blue event picks.The data have been scaled for display purposes using the rms value for each trace and filtered to 1-10 Hz.

Full-waveform event location
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.

FWM
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 6b.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 Pand S-arrivals; this ensures that the P-to-S separation is preserved.
When waveform matching is used for location and for the moment tensor inversion, we should bear in mind that they are not two separate problems.Given sufficient observations and a reasonable velocity model, we should expect that only by locating the event at the right position and by determining the right source attributes could we synthesize waveforms that match the observed ones.This is schematically illustrated in Figure 9, which compares various event focal mechanisms for which there should be only one unique solution, when location and the source mechanism are considered simultaneously.Still, our knowledge of the earth is rarely accurate enough.This means that even when we put the event at the correct location, our synthetic waveforms will not be well-aligned with our observed data, and using direct L 2 minimization for moment tensor inversion; i.e., i is the spatial derivative in the ith direction at source location s, and M ij is the ði; jÞ element of the moment tensor M.This equation becomes difficult to constrain because part of the misfit originates from velocity errors and part originates from the moment tensor solution that we are seeking to derive.Therefore, the residual discrepancy in time in each component needs to be eliminated before any L 2 -norm-based moment tensor inversion can be applied.One might think that crosscorrelation should be able to remove any nonzero residual time.However, the result from the crosscorrelation depends on the actual shape and phase of the waveforms.That is, without first having an appropriate Figure 9. Schematic representation of the full-waveform matching methodology.For a series of candidate wave trains (dashed waveforms) with known source positions and focal mechanisms, there is only one candidate event that best fits the observed data in an L 2norm sense (as indicated by the gray boxes).This match is solved by a multiattribute objective function as described in the text.moment tensor solution to generate synthetic waveforms that are like the observed data, removal of the residual time can hardly succeed using crosscorrelation-based techniques.In contrast, without proper removal of the residual time, a correct moment tensor solution is hard to obtain.To overcome this problem, we exhaustively search for all possible moment tensor solutions and for each solution, we determine its corresponding residual time for the best alignment.Then, one compares solutions from all trials and finds the optimal one, for the moment tensor solution and its related residual times.We start by applying for each candidate location r s ðx; y; zÞ an approximate origin time T init 0 correction (equation 5).Note that the origin time is the least-squares solution based on traveltime picks of the observed and synthetic waveforms.We fully take the uncertainty in phase picking into consideration; thus, we only use T init 0 as a rough estimation to align the traces for subsequent crosscorrelation-based corrections.
After the T init 0 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) 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 doublecouple focal mechanism, and we generate trial solutions for every 5°in strike, dip, and rake dimensions.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) Stack the maximum normalized crosscorrelation values for each trial solution.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 rootmean-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 T init 0 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.
The objective function for the FWM, R t , is as follows: where R w is the L 2 -waveform residual between observed and synthetic data, R m is the L 2 residual between the full moment tensor solution and many trial reference double-couple solutions, R θ is a back-azimuth constraint, and R Δτ n is the trace time shift constraint (Li et al., 2016).In this objective function, we can fully use the amplitude and phase information in the waveforms for full moment tensor inversion with the R w and R m terms.In practice, at each location ðx; y; zÞ, the first two terms R w and R m are minimized jointly by solving the linear least-squares minimization problem, with an applied stabilization factor m DC resulting in the following equation: where m ¼ ½M 11 ; M 22 ; M 33 ; M 12 ; M 13 ; M 23 T is the vector form of the moment tensor M, and I is the identity matrix.The terms R w and R m could be considered as ordinary least-squares equation (R w Þ with a generalized version of the Tikhonov regularization term R m (Phillips, 1962;Tikhonov, 1963).The terms R θ and R Δτ n are minimized separately, and then their values are combined with the first two terms to determine the best source location and mechanism.The term m DC is the moment tensor double-couple solution derived from trial comparisons of synthetic and observed waveforms.In all cases, β; γ, and η can be automatically scaled with event magnitude to match with the R w term.The term W in the R θ term is a double-cosine window, which is flatter within 1σ range of the correct minimum fit and steeper outside.Thus, it tolerates discrepancies between observed and modeled back azimuths, while penalizing back azimuths that produce large errors.The last term R Δτ n is used to penalize large crosscorrelation shifts, further improving the location accuracy.Usually R θ and R Δτ n are weak constraints, with values less than 1/10 of that of the R w term.This formulation can be generalized to accommodate alternative stabilization factors for different problems, e.g., injection or induced events.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 10 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 objec-tive 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 10.In this example, the contribution from the R M 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).

WAM versus FWM
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 11   refers to a synthetic event generated from the Green's functions and a specified moment tensor.As can be seen in Figure 11, the raybased 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.

Deep borehole network
The deep borehole network includes the SDM-1 and ZRP-1 wells (Figure 3).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 12 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 13).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 13a) and for the primary and auxiliary planes (Figure 13b).There are no obvious orientations that match the fault trends shown in Figure 12.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 14).The event location (Figure 14a) 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 14b), 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-toreceiver 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).

Shallow borehole network
The results from 100 located events using the full-waveform workflow (Appendix B) are displayed in Figure 15.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, M L , 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 15).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 16.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 17).The objective function, which is the summation of contributions from eight shallow boreholes, shows a very welldefined minimum (Figure 17a).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 17b) 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.

DISCUSSION
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 18a 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 18b.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 18c 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 18d).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 <100 m 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 19 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.

CONCLUSION
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.
Manuscript received by the Editor 19 March 2018; revised manuscript received 13 November 2018; published ahead of production 12 December 2018; published online 01 March 2019.

Figure 1 .
Figure 1.Location map of the Groningen gas field showing the distribution of earthquakes recorded since 1986 by KNMI and superimposed on the outline of the field (the gray polygon) and interpreted fault map (in blue).The events are scaled by magnitude and color coded as per the following magnitude ranges: M L < 1.5 (green), ≥ 1.5 − 2.5 (orange), and events >2.5 (red).The number of earthquakes per year for each magnitude class is shown in the graph on the right.The black line represents the total number of earthquakes recorded.Approximately two-thirds of the earthquakes detected have magnitudes that are M L < 1.5 and cannot be felt at the surface.

Figure 2 .
Figure 2. (a) Sonic log from the SDM-1 well along with the interpretation of key markers within the stratigraphic section.(b) West to east cross section from the subsurface P-wave velocity model through the location of the SDM-1 well.

Figure 3 .
Figure 3. (a) Location map of the shallow monitoring network (black dots) and deep boreholes (white triangles) superimposed on the Zechstein salt isopach.The salt is located immediately above the reservoir, and the thickness variations disrupt the earthquake energy as it travels toward the surface.(b) Deep borehole receiver configuration for the ZRP-1 and SDM-1 wells.(c) Well string geometry for a shallow borehole with geophones positioned at depths of 50, 100, 150, and 200 m below the surface.

Figure 4 .
Figure 4. Three-dimensional elastic waveform simulation demonstrating the wave propagation for three earthquake sources located at different source positions within the reservoir.The wavefield snapshots (a-c) are shown at 0.92 s from the earthquake initiation.The corresponding receiver gathers for an observer positioned at the surface are shown in (d-f) at 4 s, total time from the earthquake initiation.A double-couple source was used in each case with a maximum frequency content of 20 Hz, and the velocity model was sampled at 5 × 5 × 5 m.The dashed box in (c) indicates the location of the simulation shown in Figure 5.

Figure 5 .
Figure 5. (a) Wavefield simulation within the area defined by the box from Figure 4c, but at a higher maximum frequency of 100 Hz for a velocity model sampled at 5 × 5 × 5 m.The wavefield is shown at 0.46 s from the earthquake initiation.A double-couple source was used and is located at the left-hand boundary of the model within the reservoir.(b) Seismic displacement trace recorded at the receiver positioned within the reservoir (white inverted triangle).The most energetic arrivals observed are the P and S guided waves.
n ðtÞ is the nth component of receiver r of the observed data, G ðrÞ nj ðtÞ is the Green's function of receiver r for a vector force in the jth direction and the receiver component in the nth direction, ∂ s

Figure 8 .
Figure 8.A comparison of the key attributes (traveltime differences and polarizations) that contribute to the attribute event location objective functions.The images compare a depth slice through the objective function at a depth of 2975 m for the F 1 and F 3 terms computed for the two deep monitoring wells discussed in the main text.The summation of the F 1 and F 3 terms for each well and between wells is also shown.The summation of all the objective functions for all terms and wells is shown in the bottom right corner.The well locations are indicated by white triangles, and the blue dot indicates the final event location for this example.The color scale for each plot is shown in the top right corner and measures the normalized F 1 and F 3 attributes and their normalized sums.

Figure 10 .
Figure 10.Comparison of key terms and values from equation 8 that contribute to the full-waveform objective functions for an example event located using the two downhole wells.(a) The L 2 -waveform residual between observed and synthetic data.(b) The L 2 residual between the full moment tensor solution and many trial reference double-couple solutions.(c) Back-azimuth residual.(d) Residual time shift.(e) The summation of all terms in (a-d).The black dot indicates the final event location selected, and the black triangles are the locations of the two downhole wells, SDM-1 and ZRP-1.In each case, the parameter values have been scaled so that the variations can be observed.
Figure 11.Comparison of event locations derived using (a) raytraced traveltimes, (b) WAM workflow, and (c) FWM workflow.The red dot indicates the location derived from each workflow, and the blue dot is the true location of the event positioned on a fault.The black triangles are the locations of the two deep boreholes.The events are superimposed onto the top reservoir map, and the depth range is provided on the right.

Figure 12 .
Figure 12.Event locations for 45 events derived from the SDM-1 and ZRP-1 monitoring network (the white triangles) are shown superimposed on the top Rotliegend depth map.Note the close alignment of the events to the existing faults but also displaying a lot of variation in the focal mechanisms as indicated by the beach balls.

Figure 13 .
Figure 13.Stereonet displays for all dip and strike values calculated from the moment tensor inversion for data recorded from the downhole instrumentation, as shown in Figure 12.(a) The display on the left includes the dip and strike density contours.Black dots are the individual calculated values for each event, and the rose diagram at center shows the azimuthal variation and relative occurrence of these values.(b) The focal planes for the primary and auxiliary planes are displayed, and the rose diagram at the center shows the azimuthal variation and relative occurrence of the planes.

Figure 14
Figure 14.(a) Depth slice through the event location objective function for an event recorded by the deep borehole network.The value displayed is one minus the normalized crosscorrelation between observed and synthetic data.The wells used to constrain the solution are indicated by the black triangles, and the final location is indicated by the white dot.(b) The objective function for the focal mechanism, generated from the moment tensor inversion by scanning over a range of focal angles in dip, strike, and rake dimensions.The selected focal mechanism is indicated by the white dot.

Figure 16 .
Figure 16.Stereonet displays for all dip and strike values calculated from the moment tensor inversion for data recorded by the shallow borehole network, as shown in Figure 15.(a) The display on the left includes the density contours for the dip and strike values.Black dots are the individual calculated values for each event, and the rose diagram at center shows the azimuthal variation and relative occurrence of these values.(b) The focal planes for the primary and auxiliary planes are displayed, and the rose diagram at the center of the plot shows the azimuthal variation and relative occurrence of the focal planes.Angles are shown in degrees.

Figure 15 .
Figure 15.Event locations and focal mechanisms for 100 events recorded by the shallow borehole network between 2015 and 2017.The events are superimposed on the top Rotliegend reservoir depth map.

Figure 17 .
Figure 17.(a) Depth slice at 2950 m through the event location objective function for an event recorded by the shallow borehole network.The value displayed is one minus the normalized crosscorrelation between the observed and synthetic data.The wells used to constrain the solution are indicated by the black triangles, and the final location is shown by the white dot.(b) The objective function for the focal mechanism, generated from the moment tensor inversion by scanning over a range of focal angles in the dip, strike, and rake dimensions.The selected focal mechanism is shown by the white dot.

Figure 18 .
Figure 18.Error analysis is performed by measuring the axes of the error ellipsoid from the objective function.(a) An example depth slice through the objective function is shown for a limited azimuthal contribution from four wells (the black triangles).The displayed attribute is one minus the normalized crosscorrelation between the observed and synthetic waveforms.(b) The objective function from (a) is cropped based on the correlation value.Standard deviations are measured along the principal axes of the error ellipsoid.(c) An example depth slice through the objective function derived from 10 wells (the black triangles) that have good azimuthal coverage around the event epicenter.(d) The objective function from (c) is cropped based on the correlation value, and the standard deviations are measured along the principal axes of the error ellipsoid.

Figure 19 .
Figure 19.West to east cross section through the location of the two deep boreholes.Events located by the shallow borehole network are displayed as black dots, whereas the downhole event locations are shown in red.The figure shows the locations from the full survey area, which have been projected onto one cross section.All results have been determined to be located within the reservoir bounds.