Downloaded 09/16/17 to 24.9.229.241. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/
GEOPHYSICS, VOL. 82, NO. 6 (NOVEMBER-DECEMBER 2017); P. KS99–KS112, 14 FIGS. 10.1190/GEO2017-0263.1
Microseismic image-domain velocity inversion: Marcellus Shale case study
Ben Witten1 and Jeffrey Shragge2
image-domain adjoint-state tomography. This approach creates a complementary set of images for each chosen event through wave-equation propagation and correlating combinations of Pand S-wavefield energy. The method updates the velocity models to optimize the focal consistency of the images through adjoint-state inversion. We have determined the functionality of the method using a surface array of 192 3C geophones over a hydraulic stimulation in the Marcellus Shale. Applying the proposed joint location and velocity-inversion approach significantly improves the estimated locations. To assess the event location accuracy, we have developed a new measure of inconsistency derived from the complementary images. By this measure, the location inconsistency decreases by 75%. The method has implications for improving the reliability of microseismic interpretation with low signal-to-noise data, which may increase hydrocarbon extraction efficiency and improve risk assessment from injection-related seismicity.
ABSTRACT Seismic monitoring at injection wells relies on generating accurate location estimates of detected (micro-) seismicity. Event location estimates assist in optimizing well and stage spacings, assessing potential hazards, and establishing causation of larger events. The largest impediment to generating accurate location estimates is an accurate velocity model. For surface-based monitoring, the model should capture 3D velocity variation, yet rarely is the laterally heterogeneous nature of the velocity field captured. Another complication for surface monitoring is that the data often suffer from low signal-to-noise levels, making velocity updating with established techniques difficult due to uncertainties in the arrival picks. We use surface-monitored field data to demonstrate that a new method requiring no arrival picking can improve microseismic locations by jointly locating events and updating 3D P- and S-wave velocity models through
commercial projects (Phillips et al., 2002; Maxwell et al., 2010; Rodriguez-Pradilla, 2015) and to meet regulatory requirements (Alberta Energy Regulator, 2015; British Columbia Oil and Gas Commission, 2016). Although it is critical to have accurate location estimates, obtaining them is dependent on the accuracy of the velocity model used by the event location algorithms (Gajewski and Tessmer, 2005; Gesret et al., 2015). Ideally, the P- and S-wave velocity fields would be a full 3D model of sufficient complexity necessary to accurately locate events given the particular choice of algorithm. However, rarely does such a model exist in practice. In many circumstances, the velocity model is developed from a single or a small number of
INTRODUCTION Fluid injection into subsurface reservoirs is often monitored by a seismic network deployed to detect and locate (micro-) earthquakes. In the case of hydraulic fracturing, event locations are the primary means of evaluating the injection program, showing where new fractures have been created by the stimulation, and highlighting potential risks associated with preexisting faults (Rutledge and Phillips, 2003; Maxwell, 2014). In addition, the temporal and spatial distribution of a sequence of events can provide information about the reservoir rock and the fractures generated, such as permeability and hydraulic diffusivity (Shapiro, 2015). The information derived from microseismic monitoring is commonly used in large-scale
Manuscript received by the Editor 28 April 2017; revised manuscript received 12 June 2017; published ahead of production 21 July 2017; published online 11 September 2017. 1 The University of Western Australia, School of Earth Sciences and School of Physics, Centre for Energy Geoscience, Perth, Australia. E-mail: bwitten@ gmail.com. 2 Formerly The University of Western Australia, School of Earth Sciences and School of Physics, Centre for Energy Geoscience, Perth, Australia; presently Colorado School of Mines, Geophysics Department, Center for Wave Phenomena, Golden, Colorado, USA. E-mail:
[email protected]. © 2017 Society of Exploration Geophysicists. All rights reserved. KS99
Downloaded 09/16/17 to 24.9.229.241. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/
KS100
Witten and Shragge
well logs that may only cover the reservoir interval. Thus, there is no information about lateral heterogeneity and overburden structure (Maxwell, 2014). For borehole monitoring scenarios, the error may be limited depending on the proximity of the downhole receivers to the induced fractures. For surface monitoring, though, the deployed arrays have large areal extents, and thus a 1D assumption is often invalid. The problem becomes even more acute when there is limited velocity information above the reservoir interval and initial event location estimates are likely to be poor. Therefore, it is recommended to jointly update hypocenters and velocity estimates, rather than fix locations, to ensure optimal results (Thurber, 1992). Traditionally, earthquake hypocenter location algorithms require picking P- and S-wave arrivals. These traveltime picks are then inverted for an optimal location given a velocity model estimate (Nelson and Vidale, 1990; Moser et al., 1992). With the increased use of surface-based microseismic monitoring, new techniques based on seismic migration have been developed to locate weak events (Mw < 0) using data that generally exhibit low signal-to-noise levels (Duncan and Eisner, 2010). There are two main classes of migration event location algorithms: Kirchhoff (Kao and Shan, 2004; Baker et al., 2005) and wave equation (Artman et al., 2010; Nakata and Beroza, 2016). The results of both techniques are dependent on the input velocity model. Kirchhoff techniques calculate traveltimes by tracing rays from each point in the model space to each receiver. The data are summed across the calculated traveltime surface corresponding to each model point and a given choice of origin time. The model point with the largest stack power is then taken to be the estimated event location for the given time window. Kirchhoffbased techniques are computationally efficient, but require simplifying assumptions such as the infinite frequency approximation. Wave-equation methods use more complete wave physics to numerically back propagate a subset of the recorded data to reconstruct the wavefields. The spatial and temporal point where and when the wavefields focus is considered to be the source location and origin time. To remove the temporal dependency and avoid scanning through wavefield snapshots, it is common to apply an imaging condition to collapse the time axis. The point with the largest amplitude in the image is similarly the estimated source location of the microseismic event. This can be repeated over multiple input time windows. Although wave-equation imaging provides a more accurate representation of the reconstructed wavefields, the primary drawback is a computational expense greater than any of the other location methods. Although it is possible to provide location estimates for weak microseismic events through migration imaging, their reliability is significantly controlled by velocity model accuracy. The earthquake seismology community has developed methods of coupled hypocenter-velocity tomography (Aki and Lee, 1976; Thurber, 1983), which have more recently been applied to hydraulic fracturing (Grechka et al., 2011; Chen et al., 2017). However, one of the main drawbacks of these techniques is a reliance on picking P- and/ or S-wave arrivals in the data, which may not be possible given the weak-magnitude events commonly recorded during surface-based microseismic monitoring. Therefore, in these scenarios, it is imperative to use velocity updating procedures that do not rely on traveltime picks. One pathway to developing such a pickless velocity updating procedure is based on the knowledge that images produced from migration techniques contain velocity information. Artman et al.
(2010) note that multiple complementary zero-lag images of the same event constructed from combinations of wave modes (i.e., P-P, S-S, P-S) provide insight into the quality of the P- and S-wave velocity models. If the images do not produce consistent locations, there is likely error in the velocity model. Witten and Shragge (2015) show that in addition to the zero-lag images, one may use extended crosscorrelation images for microseismic events to investigate velocity error. Applying the self-consistency principle (i.e., all zero-lag images maximized at the same spatial location and the extended image maximum at zero lag) is a powerful qualitative means to quality control the event location estimates and produce a relative velocity update direction. Witten and Shragge (2017) extend this approach to a quantitative inversion algorithm based on the self-consistency of the suite of zero-lag and extended images. This method uses wave-equation imaging and adjoint-state tomography (Plessix, 2006) to ensure consistency among a suite of images in the zero-lag and extended image domains. Moreover, it produces P- and S-wave migration velocity model updates that provide optimized event images and location estimates. The method updates the P- and S-wave velocity model without a need for picking arrivals, making it a suitable choice for surface-based microseismic monitoring applications. Although previously validated on 2D synthetics with realistic sampling and coherent noise, this technique has yet to be implemented for microseismic field data and the associated 3D velocity model estimation. This paper aims to address this gap. In this paper, we apply the work of Witten and Shragge (2017) to a hydraulic stimulation monitoring data set from the Marcellus Shale in eastern Ohio, USA. A multiwell fracture program with an injection at approximately 1.75 km depth was monitored by a surface array consisting of 192 three-component (3C) geophones in a nominal 14 × 16 grid with spacing of approximately 400 m. Our initial information is a catalog of more than 10,000 detected events, well logs from a single vertical borehole, and basic knowledge of the geologic structure. The velocity information and geologic knowledge is representative of the extent of information commonly known prior to microseismic monitoring. We begin by summarizing the wave-equation imaging and adjoint-state inversion methodology. We then describe the field data set including the information available to us prior to starting the experiment. From these data, we show that by using a small number of detected events (28), we can update the 3D P- and S-wave velocity models to improve the image quality and the self-consistency among the suite of images. Because we do not know the true event locations, we introduce a new measure of inconsistency among the various images. Using the optimized velocity models and this measure, we subsequently image a set of 100 events to demonstrate that the velocity model improves location consistency of the events not used to constrain the inversion by 75% from 0.361 to 0.092 km. We conclude with a discussion on the quality of the resulting velocity model and event locations.
METHODOLOGY This section summarizes the joint hypocenter-velocity update procedure developed by Witten and Shragge (2017). For a more complete analysis including mathematical derivations we refer readers to this work. The image-domain adjoint-state inversion of microseismic data is an iterative five-step process: (1) wave-equation imaging, (2) image-domain residual calculation, (3) gradient calcu-
Microseismic inversion lation, (4) step-length determination, and (5) model updating. We briefly address each step in turn.
Downloaded 09/16/17 to 24.9.229.241. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/
Wave-equation imaging To produce a location estimate of a seismic event, we use a waveequation migration method. In an ideal case, we would inject recorded noise-free 3C data into the wavefield and numerically propagate it in reverse time through an accurate velocity model using the appropriate wave-equation physics to reconstruct the wavefields. The event hypocenter and origin time would be the location and time of maximum focusing in the reconstructed wavefields as determined by scanning through all wavefield snapshots (McMechan, 1982). However, in practice, we cannot exactly reconstruct the wavefields due to, for example, an inaccurate P- and/or S-wave velocity model, incomplete wavefield sampling, and incomplete physics in the numerical propagation. In addition, there is a significant computational burden of scanning through all possible origin times to locate the focus. Thus, we must examine alternate ways to determine these variables. In this study, we assume that a pseudoacoustic approximation adequately represents wave propagation, and use a wave-equation imaging technique to reduce the temporal dependence. Although an acoustic approximation does not handle the true wave modes, it is sufficient for velocity model updating and estimating event location because the overall image-domain inversion methodology is largely controlled by wavefield kinematics. In addition, the computational cost of the acoustic propagation is substantially cheaper than a full elastic scheme. Unlike an elastic propagator, using the acoustic assumption requires upfront separation of the data into P- and S-wave components rather than in the imaging condition (Witten and Shragge, 2015). We discuss data preprocessing and mode separation in the “Data description” section below. To eliminate the need to scan through wavefield propagation snapshots, we apply autocorrelation and crosscorrelation imaging conditions to stack over the time axis and produce image volumes. This reduces image-space dimensionality and improves the signal-to-noise ratio by stacking over the full source wavelet rather than at the peak amplitude. We generate three image volumes I PP , I SS , and I PS :
Z
I PP ðx; eÞ ¼
0
T
Z I SS ðx; eÞ ¼
0
T
I PS ðx; λ; eÞ ¼ max
T
0
of the P- and S-wavefields, whereas I PS is an extended image of the crosscorrelation between P- and S-wavefield energy. When λ ¼ 0 km, I PS reduces to a zero-lag image. In the following example, we only take the extended image at the location of the maximum in the zero-lag I PS image. Although it would be preferable to calculate the extended image at all spatial locations, we found a single location to be sufficient and computationally efficient when used in conjunction with the zero-lag images. We do not allow negative values in the I PS image (i.e., equation 3) to remove anticorrelated P-S energy when the wave modes are out of phase, for example, when the P-wave “catches” up and then “passes” the S-wave. This produces an artificially broad area of high-amplitude focusing, but it does not indicate colocation. In addition, we normalize each image to ensure that the residual and gradient calculations are not biased toward individual strong events. Optimal self-consistency between all three images is achieved when the maxima are colocated and the I PS image maximum is at zero lag. This scenario is taken to indicate an accurate velocity model, and the maximum is assumed to represent the correct event location. If the image maxima are inconsistent and/or an I PS maximum is not at zero lag, this indicates velocity error (Witten and Shragge, 2015), which can be used to update the velocity model. These discrepancies form the basis for the objective function described below.
Objective function and image-domain residual calculation The images generated by equations 1–3 form a set of complementary measures of the event source. Therefore, we can compare images to assess the overall quality of the velocity models. One goal is to minimize the inconsistency of the various zero-lag image pairs, which means focusing each event at a spatial location consistent among the suite of images. To measure the inconsistency in the zero-lag images, we create a penalty function Fj for each image to remove the energy around the current focal location estimate, similar to Shragge et al. (2013):
uP ðx; t; eÞuP ðx; t; eÞdt;
(1)
uS ðx; t; eÞuS ðx; t; eÞdt;
(2)
and
Z
KS101
uP ðx − λ; t; eÞuS ðx þ λ; t; eÞdt; 0 ; (3)
where T is the time length of input data; ui is the reconstructed P- or S-wavefield from the estimated P- or S-wave input data, as indicated by the subscript; e is the event index; and λ is a spatial shift between correlated wavefield components, respectively. The functions I PP and I SS are zero-lag autocorrelation images
Fj ¼ sech
αj I j ; maxðI j Þ
(4)
where αj is a parameter for the penalty width and subscript j indicates the zero-lag image type (i.e., PP, SS, and PS). This produces a penalty operator that is inversely proportional to the focal energy in the image volume. The penalty function derived from each zero-lag image can be applied to the other zero-lag images to cross-penalize and thereby remove energy where the two images are consistent. This process generates a set of image-domain residuals that can be used in the adjoint-state formalism. In the extended I PS image, we examine how well the image focuses around zero lag. This principle has been used in activesource seismology (Symes, 1993; Shen and Symes, 2008) and passive seismology (Shabelansky et al., 2015; Witten and Shragge, 2017) to invert for velocity model updates. We create a penalty function Fλ defined as a multidimensional Gaussian centered at zero lag:
Witten and Shragge
KS102
λ2y λ2 λ2 Fλ ¼ 1 − exp − x2 − − z2 ; 2σ x 2σ y 2σ z
(5)
Ki
2 X ¼ 3 Vi e
Z
0
T
∂2t ui ðx; tÞai ðx; T − tÞdtdx;
i ¼ P; S;
Downloaded 09/16/17 to 24.9.229.241. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/
(7) where σ i are the parameters used to independently adjust the width of the Gaussian function in each dimension and prevent correlated noise at far offsets from having a large impact on the residual. We omit the τ dimension because we will not apply temporal lags. After applying this penalty to the extended I PS image, the remaining energy represents the extended image-domain residual. As opposed to the more commonly used absolute value function (Shen and Symes, 2008; Yang and Sava, 2015), the Gaussian operator allows for an independent penalty to be applied in each lag direction, which is representative of the different spatial and temporal ranges over which energy focuses (e.g., due to the acquisition design). The residual energy in the penalized images is a measure of inconsistency that we use to define our objective function J over all e selected events:
J¼
1X 2 e
Z Z ϵ1
λ
−λ
F2λ ðλÞI 2PS ðλÞdλ
þ ϵ6 F2PS I 2PP þ ϵ7 F2PS I 2SS dx;
(6)
where ϵj are the weighting parameters to set the relative importance of each term. We have omitted many terms from the original objective function. Similar to Witten and Shragge (2017), we did not find they appreciably improved the result; further discussion on parameter choice can be found therein.
Gradient calculation
where ai is an adjoint wavefield. We construct the adjoint wavefields by forward propagation of an “adjoint source”, which is a function of P- and S-wavefield energy and the penalized images in equations 1–3. See Appendix A of Witten and Shragge (2017) for details.
Step-length determination and model update We first apply illumination compensation and smoothing to the calculated gradients, K P and K S . To do so we divide the gradient by the associated zero-lag image; i.e.,
K^ i ¼
Ki ; hI ii i
i ¼ P; S;
(8)
where h·i indicates multidimensional smoothing. The illumination compensation prevents the update from being concentrated at the focal point. We calculate the step length through a multiparameter line search (Tang and Ayeni, 2015). Figure 1 shows a graphical representation of the procedure. We first search for the P-wave step length, followed by that of the S-wave step length. The optimal P- and S-wave step lengths are indicated by the red and blue x’s, respectively. Using these step lengths, we perform a line search along a 2D vector whose individual components are defined by the determined individual P- and S-wave step lengths to produce a final step length. Although the gradient provides the (negative) update direction, the step length determines the appropriate magnitude of the update. We scale the gradient by the determined step length and add it to current velocity model to produce the P- and S-wave velocity model updates. All the above steps are iterated until convergence or a stopping criterion is reached.
We calculate the P- and S-wave velocity gradients, K P and K S , respectively, by means of an adjoint-state tomography formalism (see, e.g., Plessix, 2006). The negative gradients provide the direction in which to update the P- and S-wave velocity models to minimize the objective function. We obtain the gradients, K i , by computing
Geology
Figure 1. Graphical representation of a 2D step-length estimation.
The fracture job targets the Marcellus Shale Formation, which extends across much of the Appalachian Basin of the northeast United States with thicknesses between 10 and 100 m. The basin encompasses portions of numerous states including New York, Pennsylvania, West Virginia, and Ohio (Milici and Swezey, 2006). The Marcellus Shale is an organic-rich section that forms the lowest unit of the Hamilton Group, which was deposited during the middle Devonian period. The Hamilton Group is bounded by the Onondaga Limestone directly below the Marcellus Shale and the overlying Tully Limestone. The Marcellus Shale is defined by low porosity (0.5%–5.0%) and permeability (micro to nanoDarcies) (Myers, 2008), thus necessitating hydraulic fracturing to increase gas recovery and make extraction economical. The Onodaga and Tully Limestone formations act as a barrier to hydraulically induced fractures (Lee et al., 2011). Therefore, there is often limited microseismicity that is located below the Marcellus Formation (Curry et al., 2010; Ejofodomi et al., 2011; Fisher and Warpinski, 2012).
DATA DESCRIPTION
Microseismic inversion
6
X 1 (km) 4
6
Depth (km) 1 00
2
3.2 4.6 Velocity (km/s)
0
2
4
0
2
b)
4 X 1 (km)
6
6
X 2 (km)
1.8
3.2 4.6 Velocity (km/s)
6
Depth (km) 2 1 00
The initial velocity information is taken from a single dipole sonic log acquired at the well head location. The P- and S-wave velocities are measured from below the reservoir depth at approximately 1.5 km depth below mean sea level almost up to the surface. The well-log data were smoothed and extrapolated into a 3D volume, accounting for minor known structural dip in the area, approximately 2%. The velocity models have grid cells of 18 m in all directions and a total size of ðX1; X2; ZÞ ¼ ð6.42; 5.83; 3.0Þ km. Figure 3a–3c shows the initial V P and V S models and the V P ∕V S ratio, respectively. The low V P ∕V S region is the reservoir interval. Each face of the flattened cube shows a slice extracted
1.8
X 1 (km) 4
Velocity information
a)
2
We apply the method described above on a surface microseismic data set acquired during a hydraulic fracture program in eastern Ohio, USA. The hydraulic stimulation consisted of multiple wells, and more than 100 stages were completed that targeted the Marcellus Shale Formation at approximately 1.75 km below the surface. The passive seismic network design was a grid of continuously recording nodal sensors that acquired data before, during, and after the stimulation for approximately two months total recording time. The array consisted of 192 3C 4.5 Hz Sunfull geophones and Iseis Sigma digitizers distributed in a semiregular 14 × 16 rounded grid with a nominal inline and crossline spacing of 390 m, for a total area of approximately 6.5 × 6.0 km2. There is moderate elevation relief on the site, with the geophone elevations ranging between 230 and 415 m above mean sea level. Figure 2 shows the geophone locations as red dots. The dashed white box (1.5 × 1.25 km) indicates the approximate area of the stimulated volume for reference.
through the 3D volume in the X1 − X2 plane (top face), X2 − Z plane (front face), and X1 − Z plane (side face). The crosshairs on the panels indicate the extraction locations for each face. In this case, the faces shown are at Z ¼ 1.53 km, X1 ¼ 3.22 km, and X2 ¼ 2.93 km. For reference, we project approximate boundaries of the stimulated volume on the 2D faces as dashed white boxes.
0
2
4
0
2
X 1 (km) 4
c)
4 X 1 (km)
6
6
X 2 (km)
1.6
1.8
2
2.2
2
VP /VS
Depth (km) 2 1 00
Downloaded 09/16/17 to 24.9.229.241. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/
Stimulation and monitoring
KS103
0
2
4
X 2 (km)
Figure 2. Satellite image showing topography with geophone locations in red. The white box (1.5 × 1.25 km) indicates the approximate stimulated volume.
0
2
4 X 1 (km)
6
Figure 3. Initial velocity models and estimated event locations: (a) V P and I PP locations; (b) V S and I SS locations; and (c) V P ∕V S ratio and I PS locations. The dashed white boxes indicate the approximate stimulated volume.
Witten and Shragge
KS104
To improve the signal-to-noise levels, we use individual events, rather than continuous records, as input to the method described above; therefore, a detection step is necessary prior to imaging. In this case, we have access to the service company’s catalog containing more than 10,000 detected events. Of these events, we select a subset of 28 larger events that were detected during pumping from different wells and stages. We use no more than one event from any given stage to help maximize the likely spatial distribution of the events. We found that this number of events was sufficient
50
Trace 100 150
b)
Trace 100 150
Trace 100 150
1
e) 50
Trace 100 150
1
h)
i) 0 0
Trace 100 150
1
50
0
0
1
Trace 100 150
0
50
1
0
0
Easting
Vertical
g)
f)
2
2
2
1
1
0
0
0
50
0
0
0
2
2
d)
c)
Easting
Vertical
to produce reliable velocity updates and would run in a reasonable time on our compute resources. Of the selected events, the maximum, minimum, and median reported moment magnitudes are Mw ¼ 0.24, Mw ¼ −0.74, and Mw ¼ −0.45, respectively. Figure 4a–4i shows the vertical, easting, and northing components for the largest (Figure 4a–4c), smallest (Figure 4d–4f), and median (Figure 4g–4i) events after applying a band-pass filter between 5 and 80 Hz. Although the events are weak, they are discernible in the trace data. However, picking accurate individual arrivals for the P- and S-wavefield data would be difficult and prone to error, even for the largest event. We apply minimal preprocessing to the raw traces to condition the input data. After applying a 5–80 Hz band-pass filter, we take the Trace waveform envelope (Taner et al., 1979) of each 50 100 150 component. The envelope calculation primarily serves to remove polarity changes across the array due to source radiation pattern; however, it also has an advantage of lowering the frequency content of the data allowing for sparser spatial sampling of the model and larger temporal sampling during propagation. After the envelope processing, we apply a second band-pass filter between 4 and 20 Hz to return each trace of the data set to zero mean. Finally, we normalize each trace. As discussed above, we must first separate Northing the data into the P- and S-wave components for input prior to imaging. We have found it sufTrace ficient to take the P-wave data as the processed 50 100 150 vertical component, whereas the S-wave data are the average of the two processed horizontal components. Other more accurate wavefield separation methods are available (Dankbaar, 1985) and may provide modest improvements when signal-to-noise levels are high. However, no method exists to separate wave modes from very low signal-to-noise data directly in the data domain. Finally, we apply a broad windowing functioning to highlight the P- and S-wave arrivals, respectively, and we scale each trace. Figure 5a– Northing 5f shows the preprocessed P- and S-wavefield Trace data for the maximum (Figure 5a and 5b), mini50 100 150 mum (Figure 5c and 5d), and median (Figure 5e and 5f) magnitude events, respectively. These panels are the processed data from the raw traces presented in Figure 4 that are input into the inversion process.
2
2
3D INVERSION RESULTS
Vertical
2
Time (s)
50
0
0
1
1 2
Time (s)
0
0
0
a)
Time (s)
Downloaded 09/16/17 to 24.9.229.241. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/
Detected events and preprocessing
Easting
Northing
Figure 4. Vertical, easting, and northing components (left, center, and right columns) of the trace data for the (a-c) largest, (d-f) smallest, and (g-i) median magnitude events used in the inversion.
Using the initial migration velocity model estimate, we image the 28 events selected for the inversion process. To account for topographic differences in receiver location, we extend the model to the maximum receiver elevation and inject the data for each geophone at a depth relative to its elevation. For each event, we generate the I PP , I SS , and I PS images (equations 1–3), and we
X j;k¼PP;SS;PS
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m m m ðxm j − xk Þ · ðxj − xk Þ; j ≠ k;
0
Trace (s) 50 100 150
0
S Trace (s) 50 100 150
0
Trace (s) 50 100 150
0 2
1
0 1 2
Time (s)
Trace (s) 50 100 150
0 1
0
d)
P
(10)
e) 0
c0 ¼
c)
2
where c0 is the distance between the estimated event locations and ce is the distance from zero lag to the maximum in the I PS . We define c0 as
0
P Trace (s) 50 100 150
0
b)
2
(9)
Trace (s) 50 100 150
m
S
f) 0
c0 þ ce ; 4
0
a)
1
C¼
KS105
0.517 km is not shown. Event consistency is not correlated with the event magnitude at the levels we tested. Therefore, we assume that the location inconsistency is a function of velocity error rather than signal-to-noise levels. The parameters we use for the inversion are ϵ1 ¼ ϵ6 ¼ ϵ7 ¼ 1, with αPP ¼ αSS ¼ 2 for zero-lag penalty functions Fj and σ x ¼ σ y ¼ 0.180 km and σ z ¼ 0.270 km for the extended penalty function Fλ . We extend the image by ðjλx j; jλy jÞ ≤ 0.360 km in X1 and X2 and by jλz j ≤ 0.540 km in depth, with uniform
Time (s)
select the maximum in each image deeper than 0.5 km as a source location estimate. We exclude the near surface due to effects of noise and incomplete spatial sampling that combine to produce high amplitudes in the near-surface image. The white dots in Figure 3a– 3c show the event locations projected on the 2D planes for the zerolag I PP , I SS , and I PS images, respectively. We observe a wide variation of possible event locations depending on the particular choice of imaging condition. In the I PP and I SS images, there are bimodal distributions of event locations. The I PP location estimates are approximately 1.5 and 2.0 km in depth, whereas the I SS locates the events approximately 1.5 and 2.3 km in depth. The I PS -estimated event locations are more concentrated around the well depth, but they are still scattered in depth between 1.2 and 2.0 km. As mentioned above, we do not expect events to occur substantially below the reservoir interval due to the limestone formation that underlies the Marcellus Shale. Figure 6 shows images from a representative event using the initial V P and V S models. The input data are from the moderate magnitude event shown in Figure 5e and 5f. Figure 6a and 6b shows the I PP and I SS images. Figure 6c and 6d shows the zero-lag I PS image and a slice through the I PS extended volume extracted at the maximum location in the zero-lag volume. The crosshairs in all panels are at the images’ maximum location. Comparing the images, we note relatively strong focusing compared with the background energy, but inconsistency between the various source location estimates. To measure the self-consistency of the events, we propose an inconsistency value C that is the average of the distance between the estimated event location in each zero-lag image and the distance from the zero lag of the maximum in the extended image for a given event. This can be expressed as
where m is the lag location where I PS is the maxima. If all images for each event focus at the same spatial location and at zero lag in the extended imaging, the inconsistency will be 0 km. For the 28 events used in the inversion, the average initial inconsistency is C ¼ 0.104 km. Figure 7 shows the individual consistency values for each event plotted against the reported magnitude. The green circles are the values from using the initial V P and V S models (Figure 3). Note that one outlier event with an inconsistency value of
1
(11)
2
pffiffiffiffiffiffiffiffiffiffiffiffi m · m;
2
ce ¼
1
where x is the estimated source location for the image indicated by the subscript, while ce is Time (s)
Downloaded 09/16/17 to 24.9.229.241. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/
Microseismic inversion
P
S
Figure 5. Processed P- and S-wavefield data (left and right columns) for the (a and b) maximum, (c and d) minimum, and (e and f) median magnitude events used in the inversion procedure.
Witten and Shragge
lag spacing of Δλi ¼ 18 m. Figure 8a and 8b shows the percentage change from the initial velocity of the final inverted V P and V S models. Figure 9a–9c shows the inverted V P and V S models and the resulting V P ∕V S ratio, respectively. Again, we have projected the approximate well area as dashed white boxes for reference and the resulting 28 estimated event locations extracted from the I PP , I SS , and I PS images. We see from the estimated event locations that the self-consistency among the suite of images has greatly increased with most of the events falling near the well locations for all imaging conditions. Figure 10a–10d shows the zero-lag I PP , I SS , and I PS images and a slice through the I PS extended volume extracted at the maximum location in the zero-lag volume using the moderate event data (Figure 5e and 5f) and the inverted velocity model. The crosshairs intersect at the maximum location in each image. These are the same data used in Figure 6, but they are formed with the inverted rather than the initial velocity model. Comparing Figures 6 and 10, we note that the images constructed using the
3.0907
inverted model are much more consistent and the focus is stronger relative to the background energy than those presented in Figure 6. In addition, the extended image in Figure 10d has a maximum at zero lag, and it is more symmetric about this point than observed in the initial extended image (Figure 6d). Therefore, the inversion procedure and resulting velocity model has enforced self-consistency for this event. The average inconsistency using the inverted velocity model reduces to 0.049 km (less than three grid cells), which is a 53% increase of consistency compared with the estimated event locations using the initial velocity model and nearing the limit of what is resolvable by the data and image-domain tomography formalism. Figure 7 shows a comparison of the initial consistency for each event, shown as o’s and the final consistency indicated by #’s, against magnitude. A movement from right to left indicates more consistent location estimates. It is clear that by using the inverted velocity model, the inconsistency has been significantly reduced.
3.1821
3.237
1
1
3.2918
X 1 (km) 3 5
b)
X 1 (km) 3 5
a)
5
1
3 X 1 (km)
5
1
3.237
1
3 X 1 (km)
5
X1
+0.054864
(km) 0
3.3467
0.018288
1
–0.5
3.3467
X 1 (km) 3 5
5
–0.018288
d) (km) 0
c)
3 X 2 (km)
0.054864
3 X 2 (km)
2.1879
1.5843 1
Depth (km) 2 1.5 1
3.237
Depth (km) 2 1.5 1
3.2918
1
3 X 2 (km)
5
1
3 X 1 (km)
0.5
1.8038
z
Depth (km) 2 1.5 1
Downloaded 09/16/17 to 24.9.229.241. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/
KS106
5 X2
0 (km)
X1
0 (km)
Figure 6. Image volumes using the initial velocity model data of the moderate event shown in Figure 5d-5f. Panels (a-c) are the zero-lag I PP , I SS , and I PS images. Panel (d) shows a slice through the extended I PS image volume extracted at the maximum location of the zero-lag image.
Microseismic inversion
VALIDATION
–0.2
0.02
o o o # o o o oo o o o # o o o # o oo ooo o o
0.04 0.06 0.08 0.1 Inconsistency (km)
0.12
a)
0.14
Figure 7. Crossplot of inconsistency values versus calculated magnitudes Mw for the 28 events used in the inversion. Green circles show the initial inconsistency and blue #’s indicate the consistency values after the inversion. Increased consistency would move a green circle to the left.
X 1 (km) 4
0
# # #o# o # ### # # # ## o ## ## # ##
1.8
3.2 4.6 Velocity (km/s)
2
–0.8
–0.6
#
0
2
4
0
2
4 X 1 (km)
6
4
0
2
4 X 1 (km)
6
4
0
2
4 X 1 (km)
6
2
X 1 (km) 4
X 1 (km) 4
6
b)
–0.8 0 0.8 dV (percent)
1.6
Depth (km) 2 1 00
Depth (km) 2 1 0 0
2
–1.6
0
2
2
4
0
2
4 X 1 (km)
6
c)
2
X 1 (km) 4
X 1 (km) 4
6
X 2 (km)
6
X 2 (km) 0
b)
6
X 2 (km)
a)
6
Depth (km) 2 1 0 0
–0.4
Mw
o
6
We next validate the velocity model by selecting 100 random events from the catalog that were not used as constraints in the inversion. The selected events range in magnitude from −1.15 < Mw < −0.18. The events are preprocessed automatically in the same manner as the events used in the inversion.
0
0.2
#o
#
–10 0 10 dV (percent)
20
Depth (km) 2 1 0 0
2
–20
Depth (km) 2 1 00
Downloaded 09/16/17 to 24.9.229.241. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/
KS107
0
2 X 2 (km)
0
2
4 X 2 (km)
0
2
4 X 1 (km)
6
Figure 8. Percentage change in the (a) V P and (b) V S models from the 3D inversion process.
Figure 9. Inverted velocity models and 28 estimated calibrated event locations: (a) V P and I PP locations, (b) V S and I SS locations, and (c) V P ∕V S ratio and I PS locations. The dashed white boxes indicate the approximate stimulated volume.
Witten and Shragge interval (Z ¼ 1.5 km). The locations from the I SS images tend to be slightly above the reservoir, the I PP locations are below, whereas the I PS fall largely within the reservoir. Events generally do not locate significantly below the reservoir, which is unlikely due to the underlying limestone. There is a residual discrepancy within individual outlier events locating above the reservoir (I SS and I PS ), below the reservoir (I PP ), and/or laterally shifted from the wells (I PP and I PS ). This is likely due to a combination of imperfect velocity models and low signal-to-noise levels for individual events, particularly the P-wave data. The mean inconsistency has decreased from 0.485 to 0.123 km when switching from the initial to the inverted velocity model, with the median reducing from 0.361 to 0.092 km. This represents an inconsistency reduction of 75%. Figure 13a shows a crossplot of inconsistency values using the initial and inverted velocity model for the 100 validation events. The dashed line indicates no change in inconsistency. Points falling below the dashed line are events in which the inconsistency is lower using
Figure 11a–11c shows the projected estimated locations for the I PP , I SS , and I PS imaging conditions along with the projected area of stimulation on the initial V P and V S models and V P ∕V S ratio, respectively. We see a similar pattern to the initial location estimates of the 28 events used during inversion. The I PP location estimates have a bimodal distribution locating just below the reservoir interval or significantly deeper at 2.0 km. The I SS location estimates are very poor, mostly located at the reservoir, at the bottom of the model, or at the minimum allowable location of 0.5 km depth. The location estimates from the I PS images spread from the reservoir depth to 2.0 km. It is evident from Figure 11 that the event locations are very inconsistent. The mean and median inconsistency values are 0.485 and 0.361 km, respectively. Figure 12a–12c shows the I PP -, I SS -, and I PS -estimated event locations for the 100 validation events overlain on the V P and V S models and the resulting V P ∕V S ratio from the inversion process, respectively. The estimated event locations are now much more consistent in depth with the majority located around the reservoir
3.2004
3.1455
3.2187
1
1
3.2553
X 1 (km) 3 5
b)
X 1 (km) 3 5
a)
3 X 2 (km)
5
1
3 X 1 (km)
5
1
3.1455
3 X 2 (km)
1
3 X 1 (km)
5
0
X1
0
1
–0.5
3.237
X 1 (km) 3 5
5
0
d) (km) 0
c)
1.438
1.5661 1
Depth (km) 2 1.5 1
3.2187
Depth (km) 2 1.5 1
3.2553
z
0.5
1.4929
0
(km) 0
3.237 Depth (km) 2 1.5 1
Downloaded 09/16/17 to 24.9.229.241. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/
KS108
1
3 X 2 (km)
5
1
3 X 1 (km)
5 X2
0 (km)
X1
0 (km)
Figure 10. Zero-lag image volumes (a) I PP , (b) I SS , and (c) I PS using the inverted velocity model and data of the moderate event shown in Figure 5d–5f and panel (d) shows a slice through the extended I PS image volume extracted at the maximum location of the zero-lag image.
KS109
Examining the velocity updates shown in Figure 8, we note a few features. First, we stress that the resulting velocity model is an “imaging velocity” and likely does not well represent true geology. This is partly due to a limited source distribution that restricts
a)
a)
6
6
3.2 4.6 Velocity (km/s)
6
4
0
2
4 X 1 (km)
0
6
0
2
1.8
4 X 1 (km)
6
3.2 4.6 Velocity (km/s)
6
Depth (km) 2 1 00
Depth (km) 2 1 00
2
6
X 1 (km) 4
6
X 1 (km) 4
3.2 4.6 Velocity (km/s)
2
4 X 2 (km)
b)
1.8
2
6
2 X 2 (km)
2
4
0
2
X 2 (km)
4 X 1 (km)
0
6
2
2.2
1.6
4 X 1 (km)
1.8
6
2
2.2
VP /VS
2
2
VP /VS
Depth (km) 2 1 00
2
0
X 1 (km) 4
6
X 1 (km) 4
1.8
4 X 2 (km)
c)
1.6
2
6
0
c)
1.8
Depth (km) 2 1 00
Depth (km) 2 1 00 0
b)
X 1 (km) 4
3.2 4.6 Velocity (km/s)
2
1.8
DISCUSSION
2
X 1 (km) 4
6
the inverted model than the initial model. Figure 13b presents the same data, but only examines events with an inconsistency of less than 0.5 km. Although there are four events that have slightly increased inconsistency values, most have substantially decreased. This indicates that the inverted 3D velocity model improves the imaging and provides more reliable event location estimates.
Depth (km) 2 1 00
Downloaded 09/16/17 to 24.9.229.241. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/
Microseismic inversion
0
2
4 X 2 (km)
0
2
4 X 1 (km)
6
Figure 11. Initial velocity models and 100 estimated event locations: (a) V P and I PP locations; (b) V S and I SS locations; and (c) V P ∕V S ratio and I PS locations. The dashed white boxes indicate the approximate stimulated volume.
0
2
4
X 2 (km)
0
2
4 X 1 (km)
6
Figure 12. Inverted velocity models and 100 estimated event locations: (a) V P and I PP locations, (b) V S and I SS locations, and (c) V P ∕V S ratio and I PS locations. The dashed white boxes indicate the approximate stimulated volume.
Witten and Shragge
Although the event locations are much more consistent using the inverted velocity model, they are not perfectly consistent. The event location estimates from the I SS images tend to locate slightly above the wells, the I PP slightly below, and the I PS around the well depth. Therefore, there is residual error in the velocity. However, this is the optimal velocity as defined by our objective function. We may include additional terms of Witten and Shragge (2017) to reduce the inconsistency between the I PP and I SS images, but this will also most likely adversely affect the extended image. We contend that once the zero-lag images are sufficiently selfconsistent, the I PS provides the most reliable estimate of the event location. This is partly due to the inversion optimizing the I PS image through the first term. In addition, this image has higher resolution due to the differing velocities of the P- and S-wavefields and the improved signal-to-noise ratio obtained by effectively using more data to construct the image. As observed in the synthetic tests of Witten and Shragge (2017), despite having residual errors in the I PP and I SS event locations, the I PS estimates are very accurate. The inconsistency value may be used as a quality-control parameter for event location confidence. Thus, interpretations can be made on a subset of events that meet the required criteria. As an example, we set a threshold at C ¼ 0.150 km and remove all events with greater inconsistency values. Figure 14 shows the same result as Figure 12, but only for the 84 events with C ≤ 0.150 km. We see that all of the events fall near the well depths. Because all these events are self-consistent, we have greater confidence that they are accurately located. Events falling outside this confidence range should be investigated further to determine which if any of the estimated event locations may be accurate because individual components such as noisy data on the vertical component or errors in automatic windowing could cause a poor location estimate. The primary drawback of the method is its computational expense. Each imaging step requires four acoustic propagations. We use one P- and one S-wave propagation to construct the zero-lag images and another for the extended image so that we do not need to write out 4D volumes for each event. The gradient calculation requires four acoustic propagations (i.e., the forward and adjoint for the P- and S-wavefields). The 2D step length determination requires six imaging steps per event. Therefore, for a single iteration, we must calculate 18 3D acoustic wavefield propagations. Despite using graphics processing unit (GPU) propagators (Weiss and Shragge, a) b) 2013) and using an embarrassingly parallel strategy over events, each iteration took approximately a day on a single quad-GPU node. With sufficient compute resources (i.e., approximately 20 GPU nodes), the iteration cost could be brought down to the order of an hour. After velocity model determination, location estimates for individual events can be computed in three minutes. Although the computation costs may preclude this technique from currently being routinely applied, increased compute resources and ongoing research to reduce 0 1 2 3 0 0.1 0.2 0.3 0.4 0.5 the computation expense, particularly in the Initial C (km) Initial C (km) step-length calculation, may allow for application of the method in a more reasonable time Figure 13. (a) Crossplot of consistency values for 100 events using the initial and inverted velocity models. (b) A magnification only around C ≤ 0.5 km. frame. 0.5
Final C (km) 0.2 0.3 0.4 0
0.1
3
Final C (km) 1 2
the source-receiver wavepaths and, thus, limits the crossing paths that are useful for resolving structure in tomographic reconstruction. The limited source distribution also restricts the volume of the model that can be updated to improve the imaging. Therefore, events that occur far from the ones used in the inversion may still have significant location errors. In addition, using “false positive” detected events could hinder the inversion results. It is for this reason that we advocate using a subset of larger detected events that occurred during pump times at various stages and wells to increase the spatial variability and decrease the likelihood of false positives. Although the velocity field may change throughout the fluid injection due to induced fractures and stress perturbations, we have no sensitivity to these time-varying changes. This is due to using events that occur throughout the injection time frame, and thus any changes are “averaged” out. Any induced velocity perturbations would be limited to a very short portion of the entire wavepath, and thus will be difficult to resolve with the proposed methodology. We also do not observe any change in the event location accuracy, as measured by the inconsistency, with time for either the initial or final models. The inversion is solely designed to improve the focusing of the events; no geologic constraints are included through, e.g., regularization. With these caveats, we note a very large increase in the S-wave velocity at and below the reservoir depth. This increase is not geologically plausible because it decreases the V P ∕V S to 1.29, which corresponds to a negative Poisson’s ratio. The increase in the Swave velocity here is primarily due to the initial depth of the I SS locations during the first few iterations. If we assume that the final event locations are accurate, velocity updates below the event locations have little or no impact on the final locations. Therefore, this may be artifact due to a poor starting model. Another feature to note on the velocity updates is that wavepaths between the events and the far offsets are generally increased. This may be related to anisotropy in the shale layers (Sone and Zoback, 2013), which is commonly observed to affect microseismic event locations (Warpinski et al., 2009; Eisner et al., 2011). For surface data, the effect of anisotropy would be more pronounced at far offsets as the wave travels more horizontally during its cumulative path. Because the inversion uses an isotropic parameterization, we may be observing compensation to approximate the anisotropy.
0
Downloaded 09/16/17 to 24.9.229.241. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/
KS110
X 1 (km) 4
3.2 4.6 Velocity (km/s)
6
Depth (km) 2 1 00
2
1.8
0
2
4
0
2
X 1 (km) 4
b)
4 X 1 (km)
6
3.2 4.6 Velocity (km/s)
6
Depth (km) 2 1 00
2
1.8
invert the P- and S-wave velocity models using an adjoint-state method that requires no picking of arrivals making it suitable for the low signal-to-noise surface microseismic monitoring data example shown. We then validate the inverted velocity updates on a set of 100 events not used in the inversion. The velocity updates decrease the inconsistency of these events 75% from the initial model based on a proposed measure of inconsistency. The final event locations are within and above the reservoir, which is consistent with other results in the Marcellus Shale. The inverted isotropic velocity model appears to compensate for anisotropic effects at far offsets. Incorporating vertical transverse isotropy (VTI) anisotropy into the inversion is an area of future research.
We thank Statoil USA for allowing publication of the data set and ESG/Spectraseis for helping us obtain the data. We thank D. Rocha, G. Schuster, and an anonymous reviewer for their comments. B. Witten is funded by an Australian Government Research Training Program Scholarship and an Australian Society of Exploration Research Foundation grant. We acknowledge the sponsors of the UWA:RM consortium. J. Shragge acknowledges the support of Woodside Petroleum Ltd. The reproducible numeric examples use the Madagascar open-source package (http://www .reproducibility.org).
REFERENCES 0
2
4
0
2
4 X 1 (km)
6
X 1 (km) 4
6
X 2 (km)
c)
KS111
ACKNOWLEDGMENTS
6
X 2 (km)
1.6
1.8
2
2.2
2
VP /VS
Depth (km) 2 1 00
Downloaded 09/16/17 to 24.9.229.241. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/
a)
6
Microseismic inversion
0
2
4
0
X 2 (km)
2
4 X 1 (km)
6
Figure 14. Final velocity models and events with confidence value C ≤ 0.15 km: (a) V P and I PP locations, (b) V S and I SS locations, and (c) V P ∕V S ratio and I PS locations. The dashed white boxes indicate the approximate stimulated volume.
CONCLUSION We have successfully applied a recently proposed technique for updating P- and S-wave velocity models from a field microseismic data set. We use a surface-monitoring array consisting of 192 3C sensors recorded over a hydraulic fracturing job in the Marcellus Shale. From the event catalog, we use 28 microseismic events to
Aki, K., and W. H. K. Lee, 1976, Determination of three dimensional velocity anomalies under a seismic array using P-arrival times from local earthquakes. Part I: A homogeneous initial model: Journal of Geophysical Research, 81, 4381–4399, doi: 10.1029/JB081i023p04381. Alberta Energy Regulator, 2015, Subsurface Order No. 2, http://aer.ca/ documents/orders/subsurface-orders/SO2.pdf, accessed 13 March 2017. Artman, B., I. Podladtchikov, and B. Witten, 2010, Source location using time-reverse imaging: Geophysical Prospecting, 58, 861–873, doi: 10 .1111/j.1365-2478.2010.00911. Baker, T., R. Granat, and R. W. Clayton, 2005, Real-time earthquake location using Kirchhoff reconstruction: Bulletin of the Seismological Society of America, 95, 699–707, doi: 10.1785/0120040123. British Columbia Oil and Gas Commission, 2016, Ground Motion Monitoring and Submission Guideline, https://www.bcogc.ca/ground-motionmonitoring-and-submission-guideline, accessed 13 March 2017. Chen, Y., H. Zhang, Y. Miao, Y. Zhang, and Q. Liu, 2017, Back azimuth constrained double-difference seismic location and tomography for downhole microseismic monitoring: Physics of the Earth and Planetary Interiors, 264, 35–46, doi: 10.1016/j.pepi.2016.10.003. Curry, M., T. Maloney, R. A. Woodroof, and R. S. Leonard, 2010, Less sand may not be enough: Presented at the SPE Unconventional Gas Conference, SPE-131783–MS. Dankbaar, J. W. M., 1985, Separation of P- and S-waves: Geophysical Prospecting, 33, 970–986, doi: 10.1111/j.1365-2478.1985.tb00792. Duncan, P. M., and L. Eisner, 2010, Reservoir characterization using surface microseismic monitoring: Geophysics, 75, no. 5, 75A139–75A146, doi: 10.1190/1.3467760. Eisner, L., Y. Zhang, P. Duncan, M. C. Mueller, M. P. Thornton, and D. Gei, 2011, Effective VTI anisotropy for consistent monitoring of microseismic events: The Leading Edge, 30, 772–776, doi: 10.1190/1.3609092. Ejofodomi, E., J. D. Baihly, R. Malpani, R. M. Altman, T. J. Huchton, D. Welch, and J. Zieche, 2011, Integrating all available data to improve production in the Marcellus Shale: Presented at the SPE North American Unconventional Gas Conference and Exhibition, SPE-144321–MS. Fisher, M. K., and N. R. Warpinski, 2012, Hydraulic-fracture-height growth: Real data: SPE Production & Operations, 27, 8–19, doi: 10.2118/145949PA. Gajewski, D., and E. Tessmer, 2005, Reverse modelling for seismic event characterization: Geophysical Journal International, 163, 276–284, doi: 10.1111/j.1365-246X.2005.02732. Gesret, A., N. Desassis, M. Noble, T. Romary, and C. Maisons, 2015, Propagation of the velocity model uncertainties to the seismic event location: Geophysical Journal International, 200, 52–66, doi: 10.1093/gji/ggu374.
Downloaded 09/16/17 to 24.9.229.241. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/
KS112
Witten and Shragge
Grechka, V., P. Singh, and I. Das, 2011, Estimation of effective anisotropy simultaneously with locations of microseismic events: Geophysics, 76, no. 6, WC143–WC155, doi: 10.1190/geo2010-0409.1. Kao, H., and S.-J. Shan, 2004, The source-scanning algorithm: Mapping the distribution of seismic sources in time and space: Geophysical Journal International, 157, 589–594, doi: 10.1111/j.1365-246X.2004.02276. Lee, D. S., J. D. Herman, D. Elsworth, H. T. Kim, and H. S. Lee, 2011, A critical evaluation of unconventional gas recovery from the Marcellus Shale, northeastern United States: KSCE Journal of Civil Engineering, 15, 679–687, doi: 10.1007/s12205-011-0008-4. Maxwell, S., 2014, Microseismic imaging of hydraulic fracturing: Improved engineering of unconventional reservoirs: SEG. Maxwell, S., J. Rutledge, R. Jones, and M. Fehler, 2010, Petroleum reservoir characterization using downhole microseismic monitoring: Geophysics, 75, no. 5, 75A129–75A137, doi: 10.1190/1.3477966. McMechan, G. A., 1982, Determination of source parameters by wavefield extrapolation: Geophysical Journal International, 71, 613–628, doi: 10 .1111/j.1365-246X.1982.tb02788. Milici, R. C., and C. S. Swezey, 2006, Assessment of Appalachian Basin oil and gas resources: Devonian Shale-Middle and Upper Paleozoic Total Petroleum System: Technical Report: United States Geological Survey. Moser, T., T. Van Eck, and G. Nolet, 1992, Hypocenter determination in strongly heterogeneous earth models using the shortest path method: Journal of Geophysical Research, 97, 6563–6572, doi: 10.1029/91JB03176. Myers, R., 2008, Marcellus shale update: Independent Oil & Gas Association of West Virginia. Nakata, N., and G. Beroza, 2016, Reverse time migration for microseismic sources using the geometric mean as an imaging condition: Geophysics, 81, no. 2, KS51–KS60, doi: 10.1190/geo2015-0278.1. Nelson, G., and J. Vidale, 1990, Earthquake locations by 3-D finite-difference travel times: Bulletin of the Seismological Society of America, 80, 395–410. Phillips, W. S., J. T. Rutledge, L. S. House, and M. C. Fehler, 2002, Induced microearthquake patterns in hydrocarbon and geothermal reservoirs: Six case studies: Pure and Applied Geophysics, 159, 345–369, doi: 10.1007/ PL00001256. Plessix, R. E., 2006, A review of the adjoint-state method for computing the gradient of a functional with geophysical applications: Geophysical Journal International, 167, 495–503, doi: 10.1111/j.1365-246X.2006.02978. Rodriguez-Pradilla, G., 2015, Microseismic monitoring of a hydraulic-fracturing operation in a CBM reservoir: Case study in the Cerrejn Formation, Cesar-Rancheria Basin, Colombia: The Leading Edge, 34, 896–902, doi: 10.1190/tle34080896.1. Rutledge, J. T., and W. S. Phillips, 2003, Hydraulic stimulation of natural fractures as revealed by induced microearthquakes, Carthage Cotton
Valley gas field, east Texas: Geophysics, 68, 441–452, doi: 10.1190/1 .1567214. Shabelansky, A. H., A. E. Malcolm, M. C. Fehler, X. Shang, and W. L. Rodi, 2015, Source-independent full wavefield converted-phase elastic migration velocity analysis: Geophysical Journal International, 200, 952–966, doi: 10.1093/gji/ggu450. Shapiro, S., 2015, Fluid-induced seismicity: Cambridge University Press. Shen, P., and W. Symes, 2008, Automatic velocity analysis via shot profile migration: Geophysics, 73, no. 5, VE49–VE59, doi: 10.1190/1.2972021. Shragge, J., T. Yang, and P. Sava, 2013, Time-lapse image-domain tomography using adjoint-state methods: Geophysics, 78, no. 4, A29–A33, doi: 10.1190/geo2013-0044.1. Sone, H., and M. D. Zoback, 2013, Mechanical properties of shale-gas reservoir rocks. Part 1: Static and dynamic elastic properties and anisotropy: Geophysics, 78, no. 5, D381–D392, doi: 10.1190/geo2013-0050.1. Symes, W., 1993, A differential semblance criterion for inversion of multioffset seismic reflection data: Journal of Geophysical Research: Solid Earth, 98, 2061–2073, doi: 10.1029/92JB01304. Taner, M. T., F. Koehler, and R. E. Sheriff, 1979, Complex seismic trace analysis: Geophysics, 44, 1041–1063, doi: 10.1190/1.1440994. Tang, Y., and G. Ayeni, 2015, Efficient line search methods for multiparameter full wavefield inversion: U.S. Patent 20,150,323,689. Thurber, C., 1992, Hypocenter-velocity structure coupling in local earthquake tomography: Physics of the Earth and Planetary Interiors, 75, 55–62, doi: 10.1016/0031-9201(92)90117-E. Thurber, C. H., 1983, Earthquake locations and three-dimensional crustal structure in the Coyote Lake area, central California: Journal of Geophysical Research: Solid Earth, 88, 8226–8236, doi: 10.1029/ JB088iB10p08226. Warpinski, N. R., C. K. Waltman, J. Du, and Q. Ma, 2009, Anisotropy effects in microseismic monitoring: Presented at the SPE Annual Technical Conference and Exhibition, SPE-124208-MS. Weiss, R., and J. Shragge, 2013, Solving 3D anisotropic elastic wave equations on parallel GPU devices: Geophysics, 78, no. 2, F7–F15, doi: 10 .1190/geo2012-0063.1. Witten, B., and J. Shragge, 2015, Extended imaging conditions for passive seismic data: Geophysics, 80, no. 6, WC61–WC72, doi: 10.1190/ geo2015-0046.1. Witten, B., and J. Shragge, 2017, Image-domain velocity inversion and event location for microseismic monitoring: Geophysics, 82, no. 5, KS71– KS83, doi: 10.1190/geo2016-0561.1. Yang, T., and P. Sava, 2015, Image-domain wavefield tomography with extended common-image-point gathers: Geophysical Prospecting, 63, 1086–1096, doi: 10.1111/1365-2478.12204.
