ABSTRACT

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 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 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 (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 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 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.

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 tiC at station i can be written as  
tiC=T(xi,yi,zi,x0,y0,z0,)+t0,
(1)
where T is the traveltime as a function of the location of the station (xi,yi,zi) and the hypocenter. This equation has four unknowns, namely, the origin location (x0,y0,z0) and time (t0). 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 ri between the observed and calculated traveltimes, which is equivalent to the difference between the observed and calculated arrival times  
ri=ti0tiC,
(2)
where ti0 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 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.

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 4a4c 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 4d4f, 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 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 WAM objective function, F(x,y,z), to be minimized is defined using two terms of Li’s original formulation, namely, F1 describing the minimization of the observed and synthetic traveltimes and F3 describing the minimization of the polarization for the observed and synthetic data. These terms are calculated as follows:  
F(x,y,z)=λF1+(1λ)F3,
(3)
where  
F1(x,y,z)=i[(tP,iobstP,isynT0)2+(tS,iobstS,isynT0)2],
(4)
in which the difference in moveout between observed and synthetic (P and S) traveltimes is minimized and T0init is calculated to align the observed and synthetic arrival times (this can be calculated per well or for all wells simultaneously; ; depicts the average):  
T0init=0.5tPobs+tSobs0.5tPsyn+tSsyn,
(5)
and the second term,  
F3(x,y,z)=i[1(eP,isyn·eP,iobs)2],
(6)
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 (eP,iobs) and synthetic data (eP,isyn) at station i. The parameter λ in equation 1 specifies the weight assigned to the traveltime difference term F1, and it can be chosen between zero and one. The choice of weight should consider the 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 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 8. 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 8 (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.

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 P- and 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 L2 minimization for moment tensor inversion; i.e.,  
rn=13[dnobs,(r)(t)MijisGnj(r)(t)]2
(7)
will be challenging. Here, dnobs,(r)(t) is the nth component of receiver r of the observed data, Gnj(r)(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, is is the spatial derivative in the ith direction at source location s, and Mij 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 L2-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 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 rs(x,y,z) an approximate origin time T0init 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 T0init as a rough estimation to align the traces for subsequent crosscorrelation-based corrections.

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.

The objective function for the FWM, Rt, is as follows:  
Rt=Rw+Rm+Rθ+RΔτn
(8)
or  
Rt=rn=13[dnobs,(r)(tT0initΔτn(r))MijisGnj(r)(t)]2+βi=13j=i3[mijmDC,ij]2+γW(θobsθsyn)2+ηrn=13[τn(r)]2,
(9)
where Rw is the L2-waveform residual between observed and synthetic data, Rm is the L2 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 Rw and Rm terms. In practice, at each location (x,y,z), the first two terms Rw and Rm are minimized jointly by solving the linear least-squares minimization problem, with an applied stabilization factor mDC resulting in the following equation:  
m(x,y,z)=(GTG+βI)1(GTd(tT0initΔτ)+βmDC),
(10)
where m=[M11,M22,M33,M12,M13,M23]T is the vector form of the moment tensor M, and I is the identity matrix. The terms Rw and Rm could be considered as ordinary least-squares equation (Rw) with a generalized version of the Tikhonov regularization term Rm (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 mDC 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 Rw 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 Rw 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 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 10. 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).

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 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 11, 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.

RESULTS

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-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).

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, 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 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 well-defined 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.

ACKNOWLEDGMENTS

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 AND MATERIALS AVAILABILITY

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

DERIVED DEEP BOREHOLE EVENT LOCATIONS

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.

DERIVED SHALLOW BOREHOLE EVENT LOCATIONS

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.

Freely available online through the SEG open-access option.