- Email: [email protected]
Imaging hydraulic fractures by microseismic migration for downhole monitoring system
Accepted Manuscript Imaging hydraulic fractures by microseismic migration for downhole monitoring system Ye Lin, Haijiang Zhang PII: DOI: Reference:
S0031-9201(16)30109-1 http://dx.doi.org/10.1016/j.pepi.2016.06.010 PEPI 5942
To appear in:
Physics of the Earth and Planetary Interiors
Received Date: Accepted Date:
6 June 2016 21 June 2016
Please cite this article as: Lin, Y., Zhang, H., Imaging hydraulic fractures by microseismic migration for downhole monitoring system, Physics of the Earth and Planetary Interiors (2016), doi: http://dx.doi.org/10.1016/j.pepi. 2016.06.010
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
1
Imaging hydraulic fractures by microseismic migration
2
for downhole monitoring system
3 4 5
Ye Lin and Haijiang Zhang*
6 7
Wantai Microseismic Lab of School of Earth and Space Sciences
8
University of Science and Technology of China
9
96 Jinzhai Road, Hefei, Anhui 230026, China
10 11
Corresponding e-mail: [email protected]
12 13 14
Special issue on Microseismicity
15
Physics of the Earth and Planetary Interiors
16 17 18 19 20 21 22
May 23, 2016
23 24 25 26 27 1
28
ABSTRACT
29
It has been a challenge to accurately characterize fracture zones created by hydraulic
30
fracturing from microseismic event locations. This is because generally detected events
31
are not complete due to the associated low signal to noise ratio and some fracturing
32
stages may not produce microseismic events even if fractures are well developed. As a
33
result, spatial distribution of microseismic events may not well represent fractured zones
34
by hydraulic fracturing. Here, we propose a new way to characterize the fractured zones
35
by reverse time migration (RTM) of microseismic waveforms from some events. This is
36
based on the fact that fractures filled with proppants and other fluids can act as strong
37
scatterers for seismic waves. Therefore, for multi-stage hydraulic fracturing, recorded
38
waveforms from microseismic events induced in a more recent stage may be scattered by
39
fractured zones from previous stages. Through RTM of microseismic waveforms in the
40
current stage, we can determine fractured zones created in previous stages by imaging
41
area of strong scattering. We test the feasibility of this method using synthetic models
42
with different configurations of microseismic event locations and borehole sensor
43
positions for a 2D downhole microseismic monitoring system. Synthetic tests show that
44
with a few events fractured zones can be directly imaged and thus the stimulated
45
reservoir volume (SRV) can be estimated. Compared to the conventional location-based
46
SRV estimation method, the proposed new method does not depend on the completeness
47
of detected events and only a limited number of detected and located events are necessary
48
for characterizing fracture distribution. For simplicity, the 2D model is used for 2
49
illustrating the concept of microseismic RTM for imaging the fracture zone but the
50
method can be adapted to real cases in the future.
51
Keywords: Hydraulic fracturing; Downhole microseismic monitoring; Reverse time
52
migration; Stimulated reservoir volume; Scattered waves
53
1. Introduction
54
Hydraulic fracturing is an engineering tool to efficiently recover oil/gas from
55
unconventional hydrocarbon reservoirs of low-permeability, including tight sand gas/oil
56
and shale gas/oil. For hydraulic fracturing, fractures are stimulated by pumping high
57
pressure fluids into reservoirs through the treatment well. During the fracking process,
58
fractures are created by shear or tensile failure generating observed induced microseismic
59
events. Microseismic signals can be observed by surface or downhole monitoring
60
geophones. Compared to (near) surface geophones, downhole geophones installed in a
61
borehole close to the fracking well are more often used for microseismic monitoring
62
because microseismic signals have higher signal to noise ratio (SNR). Microseismic
63
events are distributed approximately around the face and tips of fractures, thus the
64
microseismic cloud provides insight into fracture location, height, length, orientation, and
65
complexity (Maxwell et al., 2010). Generally, the main purpose of microseismic
66
monitoring is to know how fractures are created, where they are distributed, and the
67
characteristics of fractures. Using microseismic monitoring, fracking design can be
68
optimized in order to improve the success ratio of hydraulic fracturing and to avoid
69
triggering large microseismic events along pre-existing faults (Maxwell, 2013). 3
70
For microseismic monitoring, the current efforts are generally focused on detecting
71
microseismic events, determining their locations and focal mechanisms. Hydraulic
72
fracture distribution, as well as stimulated reservoir volume (SRV) is generally
73
determined by microseismic locations (Urbancic et al., 1999; Mayerhofer et al., 2010).
74
The more accurate the microseismic locations, the better the fracture distribution is
75
characterized. Generally, two categories of microseismic location methodologies exist: (1)
76
arrival-based location methods and (2) migration-based location methods. In the case of
77
high SNR (Daku et al., 2004), microseismic locations are determined by picking first
78
arrivals and then finding optimal solutions by fitting absolute arrival times or arrival
79
times differences. The other category is known as the migration-based approach by
80
treating the “bright spot” as the seismic location after stacking seismic energy on each
81
node in the space domain that seismic events may occur (Kao and Shan, 2004; Drew et
82
al., 2005; Eisner et al., 2010; Chambers et al., 2010). For downhole microseismic
83
monitoring, back azimuth information of microseismic event needs to be used to
84
constrain location in addition to first arrivals. By combining microseismic locations with
85
source mechanisms, it is helpful for interpreting fracture distribution and understanding
86
fracturing process (Gharti et al., 2011; Zhebel and Eisner, 2015).
87
However, some microseismic events may be too weak to be detected and for some
88
reservoir layers hydraulic fracturing may not necessarily induce microseismic events
89
(Reshetnikov et al., 2010; Das and Zoback, 2011; Zoback et al., 2012). As a result, the
90
fracture network depicted solely by microseismic locations is biased (Johri et al., 2013; 4
91
Sicking et al., 2013; Yang and Zoback, 2014). In fact, after hydraulic fracturing, some
92
fractures are filled with proppants and become strong scatterers for seismic waves. Thus,
93
for multi-stage fracturing, fractured areas in previous stages can act as scatterers for
94
microseismic events induced by fracturing in later stages. As a result, coda waves of
95
microseismic events in later stages may contain information related to fractured zones in
96
previous stages. In this study, we will explore the possibility of using these scattered
97
waves to characterize the fractured zones.
98 99
It is known that seismic migration is the most significant process in seismic exploration because of its ability to image the subsurface structure using various kinds of
100
wavefield. Prestack depth migration has developed from Kirchhoff migration and beam
101
migration to one-way wave-equation migration and reverse time migration (Schneider,
102
1978; Hill, 1990; Claerbout, 1971; Baysal et al., 1983). In comparison, RTM is more
103
suitable for imaging complicated structures because of its ability to handle various kinds
104
of wavefield including reflected, refracted, diffracted, turning and multiple waves.
105
In recent years, the concept of seismic migration developed in seismic exploration
106
with active sources has already been applied to passive earthquake sources to image
107
structure discontinuities in the earth at various scales. Chavarria et al. (2003) used both
108
scattered P- and S-waves of recorded waveforms from earthquakes recorded in the pilot
109
well of the San Andreas Fault Observatory at Depth (SAFOD) to image the fault zone
110
structure through Kirchhoff migration. Zhang et al. (2009) used scattered P-waves of
111
earthquake waveforms to image reflectors around the SAFOD site through a generalized 5
112
Radon transform. Reshetnikov et al. (2010a) imaged the San Andreas Fault structure by
113
applying the Fresnel volume migration method to earthquake waveforms recorded by
114
downhole sensors. Reshetnikov et al. (2010b) also applied the Fresnel volume migration
115
method to induced high-frequency microseismic waveforms from the German
116
Continental Deep Drilling program to obtain the fault structure.
117
In this study, we will test the ability of using microseismic events to image fractured
118
zones by the RTM method from scattered waves in microseismic recordings. For
119
simplicity, we assume the fractures and microseismic events are located in the same
120
two-dimensional plane as the downhole sensors, but the method can be adapted to a more
121
generalized three-dimensional case. We test the feasibility of this method using different
122
configurations of microseismic event locations and sensor locations as well as different
123
synthetic models.
124
2. Microseismic reverse time migration method
125
For downhole microseismic monitoring, a string of seismic sensors is installed in a
126
borehole for monitoring microseismic events. In this case, the configuration of seismic
127
sensors and sources is similar to the vertical seismic profile (VSP) system except the
128
sources are located around reservoirs for microseismic monitoring while they are located
129
on surface for VSP. For this reason, we can adapt the VSP RTM method (Sun et al., 2011;
130
Xiao and Leaney, 2010) to downhole microseismic monitoring. The conventional reverse
131
time migration method consists of four generalized steps: (1) stimulate the source and
132
forward propagate the source wavefield, (2) set the microseismic recordings as boundary 6
133
conditions and backward propagate the wavefield, (3) apply the imaging condition to the
134
forward-propagating wavefield and the backward-propagating wavefield and get the
135
imaging result for this microseismic event, and (4) sum over imaging results from all
136
events to get the final migration image.
137
Here, we formulate the RTM method based on a 2D acoustic wave equation (Zhang
138
and Sun, 2008; Sun et al., 2011). For the borehole array, we need to compute
139
forward-propagating wavefield with event originating at , by
140 141 142 143 144 145
− ∇ , ; = 0
, ; = , −
, !
"
(1)
In the above equation, c = cx, z is the velocity model, is the microseismic source signature, and ∇ is Laplace operator.
For the common-event recordings &' , ; , ; t by borehole array at (' , ), we
need to back propagate wavefield ) based on the following equation: *
− ∇ ) , ; = 0
) = ' , ; = &' , ; , ;
"
(2)
146
To simulate both the source wavefield and the receiver wavefield, we adopted the finite
147
difference method (FDM) to solve the above equations with perfectly matched layer
148
(PML) boundary condition.
149 150 151
The imaging condition for the conventional RTM is zero-lag cross-correlation of the source wavefield and the receiver wavefield (Etgen, 1986). The images are given by +, = ! ,, , - ., , − -
7
(3)
152 153
where is time, , and . denote the forward-propagating source wavefield and backward-propagating receiver wavefield at spatial coordinate , , and +,
154
represents the imaging result at , . The final migration image is obtained by stacking
155
the images for each event.
156
The major drawback of conventional RTM is the low-frequency noise caused by
157
cross-correlation of unwanted wavefield, which adversely affects the imaging resolution.
158
Thus, the suppression of low-frequency noise is an active research area in RTM and
159
many methods have been proposed (Baysal et al., 1984; Loewenthal et al., 1987;
160
Guittonet al., 2006; Zhang and Sun, 2008; Yoon and Marfurt, 2006). In this study, we
161
adopted the wavefield-separation method (Liu and Zhang, 2007, 2011; Xiao and Leaney,
162
2010), which applies the cross-correlation imaging condition to source wavefield and
163
receiver wavefield at different directions, as follows,
164
Ix, z = ! , ± , , - .∓ , , − - -
165
Here, the superscripts “+” and “-” denote the different directions along the horizontal or
166
vertical directions.
(4)
167 168 169
3. Synthetic test of downhole microseismic migration method For simplicity, we assume that microseismic sources, downhole sensors and
170
fractures are located in a 2D plane and accordingly design a 2D acoustic model for
171
testing the microseismic migration method introduced in the above section. Furthermore,
8
172
we assume that both source locations and velocity model are known in the numerical
173
testing of the method.
174
Fig. 1 shows the configuration of microseismic events, downhole sensors and an
175
oval-shaped fracture. The fracture is located at X=1000 m and is assumed to be induced
176
by a previous hydraulic fracturing stage. Thirty-six sensors are installed in the monitoring
177
well at X=445 m. The sensors range in depth from Z=975 m to Z=2025 m with a spacing
178
of 30 m. There are 8 microseismic events, which are assumed to be induced by a later
179
fracking stage between the fracture and the monitoring well. These microseismic events
180
are randomly distributed in the fracturing zone.
181
We used a 40 Hz Ricker wavelet as the source to model synthetic waveforms based
182
on the seismic acoustic equation. The common-event recordings are 2 minutes long. A
183
10th-order accuracy central difference stencil is applied to solve the acoustic equation
184
with the perfectly matched layer (PML) boundary condition. From the snapshot of
185
forward-propagating source wavefield at t = 0.1786 s (Fig. 2a) and the recordings (Fig.
186
2b), the scattered waves from the fractured zone can be clearly seen. For microseismic
187
migration, the background velocity model is used. The respective RTM imaging results
188
for 8 individual microseismic events are shown in Fig. 3. All images clearly reveal that
189
there exists an impedance anomaly around X=1000 m where the oval-shaped fracture is
190
located. However, the images are different for different events, indicating that the
191
contribution to imaging is different for each event because of their different locations
192
relative to the fractured zone. By stacking the normalized RTM imaging result from each 9
193
individual event, the fractured zone is better imaged due to the better data coverage (Fig.
194
4). The image of the fracture shape becomes clearer and the upper tip of the fracture is
195
clearly imaged. In comparison, the lower part of the fracture is not imaged as well as the
196
upper part. This is because there are observations from only 8 events distributed in a
197
limited aperture and the sensors are located relatively higher than the fracture. As a result,
198
sensors in the monitoring well cannot receive enough wavefield that contains information
199
related to the lower part as well as the bottom tip of the fracture.
200
There are clear artifacts in the RTM imaging result for each event due to the limited
201
spatial coverage of wavefield by using only one event (Fig. 3), and the artifacts can be
202
suppressed by using 8 events (Fig. 4). There are also some artifacts appearing to the west
203
of the monitoring well (Fig.3 and Fig. 4). This kind of artifacts has also been noted when
204
applying the RTM method to vertical seismic profile system, which is mostly caused by
205
the limited spatial coverage of the observation system and can be alleviated by the
206
increase of sensor spatial coverage (Sun et al., 2011).
207
To further test how the sensor spatial coverage affects the downhole microseismic
208
RTM imaging result, we shift the sensors both downward and upward. For the sensors
209
located from Z=1575 m to Z=2625 m, the summed RTM image becomes more accurate
210
and the shape of the fracture zone becomes clearer because more wavefields coming from
211
the lower part of fracture can be received (Fig. 5a). If the sensors are shifted upward to be
212
located in the range of Z=775 m to Z=1825 m, the fracture is more poorly resolved
213
especially for the lower part of the fracture (Fig. 5b). This indicates that with a larger 10
214
aperture of sensors, the fracture can be better characterized because wavefield from
215
different directions can be received.
216
We also test how the spacing interval of borehole sensors affect the downhole
217
microseismic RTM imaging result. For comparison, we used 12 and 24 borehole sensors
218
from Z=975 m to 2025 m with an interval of 90 m and 45 m, respectively. The fracture is
219
clearly imaged from using one event or all 8 events with 12, 24 or 36 sensors (Fig. 6). It
220
is clear, however, that the fracture is slightly better imaged with more sensors. For the
221
case of using only 12 sensors, arc-shaped spatial aliasing can adversely affect the imaging
222
resolution due to the low spatial sampling of seismic wavefield. This problem can be
223
mitigated by seismic data reconstruction, such as the compressive sensing method (Tang,
224
2010).
225
In reality, the background model is not homogeneous and for the case of shale gas
226
development the model can be approximated by a layer model and the hydraulic
227
fracturing is conducted in the target shale layer. To assess how the subsurface structures
228
affect the fracture imaging, we test two models: one by adding several layers to the
229
homogeneous background (Fig. 7a), and the other by adding a normal fault in addition to
230
layers to the homogeneous background (Fig. 7c). Except for the background velocity
231
model, all the other model setup parameters are the same as those in Fig. 1. For
232
microseismic migration in this case, we again assume the background velocity model is
233
known. The tests using two models show that even with more complex structures in the
234
background model, the fracture zone can still be clearly imaged (Figure 7b, d). At the 11
235
same time the layer interfaces below the microseismic events are also imaged, especially
236
for the first layer interface (Fig. 7b, d). Furthermore, the steep normal fault to the east of
237
the fracture zone is also clearly imaged (Fig. 7d). This indicates that for the real case of
238
hydraulic fracturing, it is possible to image both fracture zone and the shale layers by
239
using microseismic events.
240
In practice, microseismic signals may be contaminated by noise. To test how noise
241
affects the imaging result, we add Gaussian distributed noise to the recordings with
242
signal-to-noise ratio (SNR) of 10 (Fig. 8a) and 5 (Fig. 8c), respectively. The normalized
243
and stacked microseismic RTM images with all 8 microseismic events show that even
244
with noise the fracture zone can still be clearly imaged although migration images are
245
contaminated to some degree (Fig. 8b, d). This demonstrates that the
246
wavefield-separation imaging condition method used in microseismic RTM behaves well
247
in suppressing the imaging noise.
248
It is inevitable that the microseismic event is generally associated with some
249
uncertainty in its location and origin time. However, with more advanced microseismic
250
location methods such as the double-difference location method (Waldhauser and
251
Ellsworth, 2000), the microseismic event locations are very accurate and the associated
252
uncertainties could be on the order of 10 m (Castellanos and Van, 2013; Li et al., 2013).
253
Here, we further test how location errors affect the RTM imaging results by perturbing
254
event locations randomly within 10 meters (Fig. 9a), 20 meters (Fig. 9c) and 50 meters
255
(Fig. 9e), respectively. Fig. 9b, Fig. 9d and Fig. 9f show the normalized and stacked RTM 12
256
images corresponding to perturbed event locations in Fig. 9a, Fig. 9c and Fig. 9e,
257
respectively. Compared with the image using true event locations shown in Fig. 4, the
258
RTM image is almost not affected by randomly perturbing locations within 10 meters
259
(Fig. 9b). In comparison, greater location errors may adversely affect the RTM result and
260
it can be seen that the fracture structures are slightly distorted by randomly perturbing
261
locations within 20 meters (Fig. 9d). If we further perturb event locations randomly
262
within 50 m on the order of the half wavelength (Fig. 9e), the fracture zone is more
263
distorted but it is still clearly recognizable in the migrated image (Fig. 9f).
264 265 266
4. Discussion and conclusions In this study, we propose to take advantage of microseismic scattered waves and
267
apply the RTM method to directly image the fracture zone by hydraulic fracturing.
268
Unlike the conventional location-based fracture characterization method that uses
269
microseismic locations as a proxy for determining the fracture distribution, the
270
microseismic migration method that we propose to use in this study does not need to
271
detect and use all microseismic events to accurately characterize the fracture zone. As
272
long as a few microseismic events with high SNR are detected and located, they can be
273
used to characterize the fractured zone. The numerical tests showed that with only 8
274
microseismic events, we are able to clearly characterize the fracture zone.
275 276
For microseismic monitoring, it is important to estimate the SRV. More importantly, it is better to estimate the effective SRV where the proppants exist and fractures are 13
277
conductively connected (Mayerhofer et al., 2010; Johri et al., 2013). However, for the
278
zone where the microseismic events are detected and located, it does not necessarily
279
indicate the existence of proppants and opened fractures there. Therefore, it has been a
280
challenge to use microseismic monitoring to characterize where proppants migrate. With
281
the microseismic migration method, it is possible to differentiate the region with or
282
without proppants and opened fractures. Even if the region is associated with
283
microseismic events, if it does not show strong scattering effect the fractures may be
284
closed due to the lack of support by proppants.
285
We showed that for different sensor positions relative to the fracture zone, the RTM
286
imaging result is different. The monitoring system having better spatial sensor coverage
287
can better image the fracture zone. Therefore, if it is possible, the sensors should be
288
installed around the same depth zone of fracking. It is also noted that different sensor
289
spacing intervals also affect the imaging result and for the case of large spacing interval
290
the spatial aliasing caused by low spatial sampling of the wavefield will impart the
291
imaging resolution. In addition to imaging the fracture zone, from different synthetic tests
292
it is also noted that the backward-propagating wavefield is focused at the source locations
293
(Fig. 3).
294
In real cases, the microseismic signals are contaminated with noise to some degrees.
295
The synthetic tests by adding different levels of noise to microseismic recordings show
296
that the wavefield-separation imaging condition method used in microseismic RTM can
297
satisfactorily deal with the noisy data to recover the fractured zone. We also test how 14
298
microseismic event location errors affect the RTM imaging results and it shows that the
299
fracture zone can be biased to some degree but can still be clearly imaged even when the
300
event locations are randomly perturbed within 50 meters, approximately on the order of
301
half wavelength. With more advanced microseismic location methods available, the
302
location errors could well be on the order of 10s of meters. The uncertainty in the origin
303
time bears a similar effect on the RTM imaging result as location errors. It is also
304
possible that different microseismic events may have different source time functions.
305
However, due to the nature of high frequency and limited bandwidth for microseismic
306
signals, the corresponding source time functions can be well approximated by a smooth
307
ramp function with a very short duration (Li et al., 2011a, b). In addition, most of
308
microseismic events from the same fracking stage generally have similar focal
309
mechanisms. Therefore, it is expected that source signatures of microseismic events
310
would not greatly affect the RTM image.
311
In these tests, we assume the background velocity model used for microseismic
312
migration is known. In real cases, the velocity model can be well estimated from sonic
313
well logging and calibrated by perforation shots (Zhang et al., 2016). Furthermore, it can
314
also be estimated by microseismic velocity tomography using arrival times from
315
microseismic events (Zhang et al., 2009). In this study, we only use the reverse time
316
migration method to the scattered P waves, but the method can also be adapted to use
317
scattered S waves based on the elastic wave equation. We emphasize that the goal of this
318
study is to show the concept of microseismic RTM to image the fractures. For the 15
319
application of this method to the real microseismic data, we need to adapt the
320
microseismic RTM to the 3D case. One way to do this is to use the 3D acoustic equation
321
and modify the equations (1) and (2) accordingly to realize the 3D microseismic RTM.
322
Another way is to apply the 2D RTM to each microseismic event that forms a 2D plane
323
with the vertical monitoring array. By stacking RTM images from different microseismic
324
events, we may be able to obtain the 3D RTM image. To better take care of uncertainties
325
in the origin time, location and source signature of microseismic event, we could apply
326
the elastic reverse time migration imaging with both microseismic scattered P and S
327
waves (Xiao and Leaney, 2010). When combined with microseismic locations and focal
328
mechanisms, microseismic migration can help us to better understand the fracture
329
network, and thus to better optimize the fracturing design and to enhance the oil/gas
330
recovery.
331 332 333
Acknowledgements This research is supported by Natural Science Foundation of China under the grant
334
No. 41274055. We are grateful for constructive comments from two anonymous
335
reviewers that are helpful for improving this paper.
336
References
337
Baysal E, Kosloff D D, Sherwood J W C. Reverse time migration[J]. Geophysics, 1983,
338
48(11): 1514-1524. http://dx.doi.org/10.1190/1.1441434
16
339 340 341 342 343
Baysal E, Kosloff D D, Sherwood J W C. A two-way non reflecting wave equation: Geophysics, 1984, 49, 132-141. http://dx.doi.org/10.1190/1.1441644 Castellanos F, Van der Baan M. Microseismic event locations using the double-difference algorithm[J]. CSEG Recorder, 2013, 38(3): 26-37. Chambers K, Kendall J, Brandsberg-Dahl S, et al. Testing the ability of surface arrays to
344
monitor microseismic activity[J]. Geophysical Prospecting, 2010, 58(5): 821-830.
345
http://dx.doi.org/10.1111/j.1365-2478.2010.00893.x
346
Chavarria J A, Malin P, Catchings R D, et al. A look inside the San Andreas fault at
347
Parkfield through vertical seismic profiling[J]. Science, 2003, 302(5651): 1746-1748.
348
http://dx.doi.org/10.1126/science.1090711
349
Daku B L F, Salt J E, Sha L. An algorithm for locating microseismic events[C]. Electrical
350
and Computer Engineering, 2004. Canadian Conference on. IEEE, 2004, 4:
351
2311-2314. http://dx.doi.org/10.1109/CCECE.2004.1347708
352
Das I, Zoback M D. Long-period, long-duration seismic events during hydraulic fracture
353
stimulation of a shale gas reservoir[J]. SEG Technical Program Expanded Abstracts
354
2011: pp. 1473-1477. http://dx.doi.org/10.1190/1.3627480
355
Drew J, Leslie H, Armstrong P, et al. Automated Microseismic Event Detection and
356
Location by Continuous Spatial Mapping[C]. SPE Annual Technical Conference and
357
Exhibition. 2005. http://dx.doi.org/10.2118/95513-MS
358
Eisner L, Hulsey B J, Duncan P, et al. Comparison of surface and borehole locations of
359
induced seismicity[J]. Geophysical Prospecting, 2010, 58(5): 809-820.
360
http://dx.doi.org/10.1111/j.1365-2478.2010.00867.x
361
Gharti H N, Oye V, Kühn D, et al. Simultaneous microearthquake location and
362
moment-tensor estimation using time-reversal imaging[C]. 2011 SEG Annual
363
Meeting. Society of Exploration Geophysicists, 2011.
364
http://dx.doi.org/10.1190/1.3627516
17
365
Guitton A, Kaelin B, Biondi B. Least-square attenuation of reverse time migration
366
artifacts[C]. SEG Technical Program Expanded Abstracts 2006, pp. 2348-2352.
367
http://dx.doi.org/10.1190/1.2370005
368 369 370
Hill N R. Gaussian beam migration[J]. Geophysics, 1990, 55(11): 1416-1428. http://dx.doi.org/10.1190/1.1442788 Johri M, Zoback M D. The Evolution of Stimulated Reservoir Volume during Hydraulic
371
Stimulation of Shale Gas Formations. Unconventional Resources Technology
372
Conference, Denver, Colorado, 12-14 August 2013: pp. 1661-1671.
373
http://dx.doi.org/10.1190/urtec2013-170
374
Kao H, Shan S J. The source-scanning algorithm: Mapping the distribution of seismic
375
sources in time and space[J]. Geophysical Journal International, 2004, 157(2):
376
589-594. http://dx.doi.org/10.1111/j.1365-246X.2004.02276.x
377
Li J, Sadi Kuleli H, Zhang H, et al. Focal mechanism determination of induced
378
microearthquakes in an oil field using full waveforms from shallow and deep seismic
379
networks[J]. Geophysics, 2011a, 76(6): WC87-WC101. http://dx.doi.org/doi:
380
10.1190/geo2011-0030.1
381
Li J, Zhang H, Kuleli H S, et al. Focal mechanism determination using high-frequency
382
waveform matching and its application to small magnitude induced earthquakes[J].
383
Geophysical Journal International, 2011b, 184(3): 1261-1274.
384
http://dx.doi.org/10.1111/j.1365-246X.2010.04903.x
385
Li J, Zhang H, Rodi W L, et al. Joint microseismic location and anisotropic tomography
386
using differential arrival times and differential backazimuths[J]. Geophysical Journal
387
International, 2013, 195(3): 1917-1931. http://dx.doi.org/10.1093/gji/ggt358
388
Liu F, Zhang G, Morton S, et al. Reverse-time migration using one-way wavefield
389
imaging condition, 77th Annual International Meeting, SEG, 2007, pp. 2170-2174.
390
http://dx.doi.org/10.1190/1.2792917
18
391
Liu, F, Zhang G, Morton S, et al. An effective imaging condition for reverse-time
392
migration using wavefield decomposition: Geophysics, 2011, 76, s29-s39.
393
http://dx.doi.org/10.1190/1.3533914
394
Loewenthal D, Stoffa P A, Faria E L. Suppressing the unwanted reflections of the full
395
wave equation: Geophysics, 1987, 52, 1007-1012.
396
http://dx.doi.org/10.1190/1.1442352
397
Maxwell S C, Rutledge J, Jones R, et al. Petroleum reservoir characterization using
398
downhole microseismic monitoring[J]. Geophysics, 2010, 75(5): 75A129-75A137.
399
http://dx.doi.org/10.1190/1.3477966
400
Maxwell S C, Urbancic T I. The role of passive microseismic monitoring in the
401
instrumented oil field[J]. The Leading Edge, 2001, 20(6): 636-639.
402
http://dx.doi.org/10.1190/1.1439012
403 404 405
Maxwell S. Unintentional seismicity induced by hydraulic fracturing[J]. Canadian Society of Exploration Geophysics Recorder, 2013, 38(8): 40-49. Mayerhofer M, Lolon E, Warpinski N, et al. What Is Stimulated Reservoir Volume?[J].
406
SPE Production & Operations, 2010, 25(1): 89-98.
407
http://dx.doi.org/10.2118/119890-PA
408
Reshetnikov A, Buske S, Shapiro S A. Seismic imaging using microseismic events:
409
Results from the San Andreas Fault System at SAFOD[J]. Journal of Geophysical
410
Research: Solid Earth (1978–2012), 2010a, 115(B12).
411
http://dx.doi.org/10.1029/2009JB007049
412
Reshetnikov A, Kummerow J, Buske S, et al. Microseismic imaging from a single
413
geophone: KTB[C]//2010 SEG Annual Meeting. Society of Exploration
414
Geophysicists, 2010b. http://dx.doi.org/10.1190/1.3513252
415 416
Schneider W A. Integral formulation for migration in two and three dimensions[J]. Geophysics, 1978, 43(1): 49-76. http://dx.doi.org/10.1190/1.1440828
19
417
Sicking C, Vermilye J, Geiser P, et al. Permeability field imaging from microseismic[J].
418
SEG Technical Program Expanded Abstracts 2012: pp. 1-5.
419
http://dx.doi.org/10.1190/segam2012-1383.1
420
Sun W B, Sun Z D, Zhu X H. Amplitude preserved VSP reverse time migration for
421
angle-domain CIGs extraction[J]. Applied Geophysics, 2011, 8(2): 141-149.
422
http://dx.doi.org/10.1007/s11770-011-0277-1
423 424 425
Tang G. Seismic data reconstruction and denoising based on compressive sensing and sparse representation[J]. Beijing: Tsinghua University, 2010 Urbancic T I J, Shumila V J, Rutledge J T J, et al. Determining hydraulic fracture
426
behavior using microseismicity[C]. Vail Rocks 1999, The 37th US Symposium on
427
Rock Mechanics (USRMS). American Rock Mechanics Association, 1999.
428
Waldhauser F, Ellsworth W L. A double-difference earthquake location algorithm:
429
Method and application to the northern Hayward fault, California[J]. Bulletin of the
430
Seismological Society of America, 2000, 90(6): 1353-1368.
431
http://dx.doi.org/10.1785/0120000006
432
Xiao X, Leaney W S. Local vertical seismic profiling (VSP) elastic reverse-time
433
migration and migration resolution: Salt-flank imaging with transmitted P-to-S
434
waves[J]. Geophysics, 2010, 75(2): S35-S49. http://dx.doi.org/10.1190/1.3309460
435
Yang Y, Zoback M D. The role of preexisting fractures and faults during multistage
436
hydraulic fracturing in the Bakken Formation[J]. Interpretation, 2014, 2(3):
437
SG25-SG39. http://dx.doi.org/10.1190/INT-2013-0158.1
438 439 440
Yoon K, Marfurt K J. Reverse-time migration using the poynting vector [J]: Exploration Geophysics, 2006, 59, 102-107. http://dx.doi.org/10.1071/EG06102 Zhang H, Wang P, van der Hilst R D, et al. Three-dimensional passive seismic waveform
441
imaging around the SAFOD site, California, using the generalized Radon
442
transform[J]. Geophysical Research Letters, 2009, 36(23).
443
http://dx.doi.org/10.1029/2009gl040372 20
444
Zhang, H., S. Sarkar, M.N. Toksoz, H.S. Kuleli, and F. Al-Kindy, Passive seismic
445
tomography using induced seismicity at a petroleum field in Oman, Geophysics, 2009,
446
74(6), WCB67, doi:10.1190/1.3253059.
447
Zhang, Jiewen, H. Zhang, Y. Zhang, and Q. Liu, Fast one-dimensional velocity model
448
determination for microseismic monitoring using station-pair differential arrival times
449
based on the differential evolution method, Physics of the Earth and Planetary
450
Interiors, in revision 2016.
451
Zhang Y, Sun J. Practical issues of reverse time migration-true-amplitude gathers, noise
452
removal and harmonic-source encoding[C]. 70th EAGE Conference & Exhibition.
453
2008. http://dx.doi.org/10.3997/2214-4609.20147708
454
Zhebel O, Eisner L. Simultaneous microseismic event localization and source mechanism
455
determination[J]. Geophysics, 2014, 80(1): KS1-KS9.
456
http://dx.doi.org/10.1190/segam2012-1033.1
457
Zoback M D, Kohli A, Das I, et al. The importance of slow slip on faults during
458
hydraulic fracturing stimulation of shale gas reservoirs[C]. SPE Americas
459
Unconventional Resources Conference. Society of Petroleum Engineers, 2012.
460
http://dx.doi.org/10.2118/155476-MS
461 462 463 464 465 466 467 468 469 470 21
471
Figure Captions
472 473
Fig. 1. 2D P-wave velocity model with a fracture zone at X=1000 m. The background
474
velocity is 4200 m/s and the fracture zone is associated with a low velocity anomaly of 15%
475
lower than the background velocity. Monitoring well is located at X=445 m with 36
476
sensors (denoted by triangles) ranging in depth from Z=975 m to Z=2025 m with a
477
spacing of 30 m. Microseismic events are (denoted by asterisks) located at (850 m, 2120
478
m), (850 m, 1975 m), (850 m, 2260 m), (862 m, 2048 m), (862 m, 2170 m), (830 m, 2107
479
m), (843 m, 2041 m), and (842 m, 2235 m), respectively.
480 481
Fig. 2. (a) Snapshot of source wavefield at t=0.1786 s. (b) Common-shot recordings for
482
the event located at (850 m, 1975 m) and the sample interval is 0.595 ms. The model
483
setup is shown in Fig. 1.
484 485
Fig. 3. RTM imaging result from 8 individual microseismic events (denoted by asterisks)
486
for the model setup in Fig. 1.
487 488
Fig. 4. Normalized and stacked RTM imaging result from 8 microseismic events shown
489
in Fig. 3.
490 491
Fig. 5. Comparison of RTM imaging results for sensors located at different depths for the
492
model setup shown in Fig. 1. (a) From Z=1575 m to 2625 m. (b) From Z=775 m to 1825
493
m.
494 495
Fig. 6. Comparison of imaging results using (left) the individual event located at (850 m,
496
2120 m) and (right) all 8 events together in different cases of having (top) 12 sensors,
497
(middle) 24 sensors and (bottom) 36 sensors. The span of these sensors is the same in 22
498
depth but with different spacing of 90, 45 and 30 m from top to bottom. Other symbols
499
are the same as those in Fig. 3.
500 501
Fig. 7 Comparison of the RTM images when the background velocity models are
502
complex. (a) The background velocity model contains several layers with different
503
velocities compared to that shown in Fig. 1; (b) Normalized and stacked RTM image
504
from 8 microseismic events for the model setup in (a); (c) The same as (a) but with a
505
normal fault to the east of the fracture zone; (d) Normalized and stacked RTM image for
506
the model setup in (c). The model layers are denoted by dashed lines in (b) and (d).
507 508
Fig. 8 The RTM images from microseismic recordings with different levels of noise. (a)
509
The recordings from a microseismic event by adding Gaussian distributed noise with
510
SNR of 10. (b) Normalized and stacked RTM image for 8 events for the noise level in (a).
511
(c) The recordings from a microseismic event by adding Gaussian distributed noise with
512
SNR of 5. (d) Normalized and stacked RTM image for 8 events for the noise level in (c).
513
The model setup is the same as Fig. 1.
514 515
Fig. 9 Normalized and stacked RTM images (right) from 8 microseismic events when
516
their locations are randomly perturbed at different degrees (left). (a) and (b) show
517
locations randomly perturbed within 10 m and the corresponding RTM image. (c) and (d)
518
show locations randomly perturbed within 20 m and the corresponding RTM image. (e)
519
and (f) show locations randomly perturbed within 50 m and the corresponding RTM
520
image. White asterisks denote true seismic locations and red asterisks denote perturbed
521
locations.
522
23
523 524
24
525 526
25
527 528
26
529 530
27
531 532
28
533 534
29
535 536
30
537 538
31
539
32
S0031-9201(16)30109-1 http://dx.doi.org/10.1016/j.pepi.2016.06.010 PEPI 5942
To appear in:
Physics of the Earth and Planetary Interiors
Received Date: Accepted Date:
6 June 2016 21 June 2016
Please cite this article as: Lin, Y., Zhang, H., Imaging hydraulic fractures by microseismic migration for downhole monitoring system, Physics of the Earth and Planetary Interiors (2016), doi: http://dx.doi.org/10.1016/j.pepi. 2016.06.010
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
1
Imaging hydraulic fractures by microseismic migration
2
for downhole monitoring system
3 4 5
Ye Lin and Haijiang Zhang*
6 7
Wantai Microseismic Lab of School of Earth and Space Sciences
8
University of Science and Technology of China
9
96 Jinzhai Road, Hefei, Anhui 230026, China
10 11
Corresponding e-mail: [email protected]
12 13 14
Special issue on Microseismicity
15
Physics of the Earth and Planetary Interiors
16 17 18 19 20 21 22
May 23, 2016
23 24 25 26 27 1
28
ABSTRACT
29
It has been a challenge to accurately characterize fracture zones created by hydraulic
30
fracturing from microseismic event locations. This is because generally detected events
31
are not complete due to the associated low signal to noise ratio and some fracturing
32
stages may not produce microseismic events even if fractures are well developed. As a
33
result, spatial distribution of microseismic events may not well represent fractured zones
34
by hydraulic fracturing. Here, we propose a new way to characterize the fractured zones
35
by reverse time migration (RTM) of microseismic waveforms from some events. This is
36
based on the fact that fractures filled with proppants and other fluids can act as strong
37
scatterers for seismic waves. Therefore, for multi-stage hydraulic fracturing, recorded
38
waveforms from microseismic events induced in a more recent stage may be scattered by
39
fractured zones from previous stages. Through RTM of microseismic waveforms in the
40
current stage, we can determine fractured zones created in previous stages by imaging
41
area of strong scattering. We test the feasibility of this method using synthetic models
42
with different configurations of microseismic event locations and borehole sensor
43
positions for a 2D downhole microseismic monitoring system. Synthetic tests show that
44
with a few events fractured zones can be directly imaged and thus the stimulated
45
reservoir volume (SRV) can be estimated. Compared to the conventional location-based
46
SRV estimation method, the proposed new method does not depend on the completeness
47
of detected events and only a limited number of detected and located events are necessary
48
for characterizing fracture distribution. For simplicity, the 2D model is used for 2
49
illustrating the concept of microseismic RTM for imaging the fracture zone but the
50
method can be adapted to real cases in the future.
51
Keywords: Hydraulic fracturing; Downhole microseismic monitoring; Reverse time
52
migration; Stimulated reservoir volume; Scattered waves
53
1. Introduction
54
Hydraulic fracturing is an engineering tool to efficiently recover oil/gas from
55
unconventional hydrocarbon reservoirs of low-permeability, including tight sand gas/oil
56
and shale gas/oil. For hydraulic fracturing, fractures are stimulated by pumping high
57
pressure fluids into reservoirs through the treatment well. During the fracking process,
58
fractures are created by shear or tensile failure generating observed induced microseismic
59
events. Microseismic signals can be observed by surface or downhole monitoring
60
geophones. Compared to (near) surface geophones, downhole geophones installed in a
61
borehole close to the fracking well are more often used for microseismic monitoring
62
because microseismic signals have higher signal to noise ratio (SNR). Microseismic
63
events are distributed approximately around the face and tips of fractures, thus the
64
microseismic cloud provides insight into fracture location, height, length, orientation, and
65
complexity (Maxwell et al., 2010). Generally, the main purpose of microseismic
66
monitoring is to know how fractures are created, where they are distributed, and the
67
characteristics of fractures. Using microseismic monitoring, fracking design can be
68
optimized in order to improve the success ratio of hydraulic fracturing and to avoid
69
triggering large microseismic events along pre-existing faults (Maxwell, 2013). 3
70
For microseismic monitoring, the current efforts are generally focused on detecting
71
microseismic events, determining their locations and focal mechanisms. Hydraulic
72
fracture distribution, as well as stimulated reservoir volume (SRV) is generally
73
determined by microseismic locations (Urbancic et al., 1999; Mayerhofer et al., 2010).
74
The more accurate the microseismic locations, the better the fracture distribution is
75
characterized. Generally, two categories of microseismic location methodologies exist: (1)
76
arrival-based location methods and (2) migration-based location methods. In the case of
77
high SNR (Daku et al., 2004), microseismic locations are determined by picking first
78
arrivals and then finding optimal solutions by fitting absolute arrival times or arrival
79
times differences. The other category is known as the migration-based approach by
80
treating the “bright spot” as the seismic location after stacking seismic energy on each
81
node in the space domain that seismic events may occur (Kao and Shan, 2004; Drew et
82
al., 2005; Eisner et al., 2010; Chambers et al., 2010). For downhole microseismic
83
monitoring, back azimuth information of microseismic event needs to be used to
84
constrain location in addition to first arrivals. By combining microseismic locations with
85
source mechanisms, it is helpful for interpreting fracture distribution and understanding
86
fracturing process (Gharti et al., 2011; Zhebel and Eisner, 2015).
87
However, some microseismic events may be too weak to be detected and for some
88
reservoir layers hydraulic fracturing may not necessarily induce microseismic events
89
(Reshetnikov et al., 2010; Das and Zoback, 2011; Zoback et al., 2012). As a result, the
90
fracture network depicted solely by microseismic locations is biased (Johri et al., 2013; 4
91
Sicking et al., 2013; Yang and Zoback, 2014). In fact, after hydraulic fracturing, some
92
fractures are filled with proppants and become strong scatterers for seismic waves. Thus,
93
for multi-stage fracturing, fractured areas in previous stages can act as scatterers for
94
microseismic events induced by fracturing in later stages. As a result, coda waves of
95
microseismic events in later stages may contain information related to fractured zones in
96
previous stages. In this study, we will explore the possibility of using these scattered
97
waves to characterize the fractured zones.
98 99
It is known that seismic migration is the most significant process in seismic exploration because of its ability to image the subsurface structure using various kinds of
100
wavefield. Prestack depth migration has developed from Kirchhoff migration and beam
101
migration to one-way wave-equation migration and reverse time migration (Schneider,
102
1978; Hill, 1990; Claerbout, 1971; Baysal et al., 1983). In comparison, RTM is more
103
suitable for imaging complicated structures because of its ability to handle various kinds
104
of wavefield including reflected, refracted, diffracted, turning and multiple waves.
105
In recent years, the concept of seismic migration developed in seismic exploration
106
with active sources has already been applied to passive earthquake sources to image
107
structure discontinuities in the earth at various scales. Chavarria et al. (2003) used both
108
scattered P- and S-waves of recorded waveforms from earthquakes recorded in the pilot
109
well of the San Andreas Fault Observatory at Depth (SAFOD) to image the fault zone
110
structure through Kirchhoff migration. Zhang et al. (2009) used scattered P-waves of
111
earthquake waveforms to image reflectors around the SAFOD site through a generalized 5
112
Radon transform. Reshetnikov et al. (2010a) imaged the San Andreas Fault structure by
113
applying the Fresnel volume migration method to earthquake waveforms recorded by
114
downhole sensors. Reshetnikov et al. (2010b) also applied the Fresnel volume migration
115
method to induced high-frequency microseismic waveforms from the German
116
Continental Deep Drilling program to obtain the fault structure.
117
In this study, we will test the ability of using microseismic events to image fractured
118
zones by the RTM method from scattered waves in microseismic recordings. For
119
simplicity, we assume the fractures and microseismic events are located in the same
120
two-dimensional plane as the downhole sensors, but the method can be adapted to a more
121
generalized three-dimensional case. We test the feasibility of this method using different
122
configurations of microseismic event locations and sensor locations as well as different
123
synthetic models.
124
2. Microseismic reverse time migration method
125
For downhole microseismic monitoring, a string of seismic sensors is installed in a
126
borehole for monitoring microseismic events. In this case, the configuration of seismic
127
sensors and sources is similar to the vertical seismic profile (VSP) system except the
128
sources are located around reservoirs for microseismic monitoring while they are located
129
on surface for VSP. For this reason, we can adapt the VSP RTM method (Sun et al., 2011;
130
Xiao and Leaney, 2010) to downhole microseismic monitoring. The conventional reverse
131
time migration method consists of four generalized steps: (1) stimulate the source and
132
forward propagate the source wavefield, (2) set the microseismic recordings as boundary 6
133
conditions and backward propagate the wavefield, (3) apply the imaging condition to the
134
forward-propagating wavefield and the backward-propagating wavefield and get the
135
imaging result for this microseismic event, and (4) sum over imaging results from all
136
events to get the final migration image.
137
Here, we formulate the RTM method based on a 2D acoustic wave equation (Zhang
138
and Sun, 2008; Sun et al., 2011). For the borehole array, we need to compute
139
forward-propagating wavefield with event originating at , by
140 141 142 143 144 145
− ∇ , ; = 0
, ; = , −
, !
"
(1)
In the above equation, c = cx, z is the velocity model, is the microseismic source signature, and ∇ is Laplace operator.
For the common-event recordings &' , ; , ; t by borehole array at (' , ), we
need to back propagate wavefield ) based on the following equation: *
− ∇ ) , ; = 0
) = ' , ; = &' , ; , ;
"
(2)
146
To simulate both the source wavefield and the receiver wavefield, we adopted the finite
147
difference method (FDM) to solve the above equations with perfectly matched layer
148
(PML) boundary condition.
149 150 151
The imaging condition for the conventional RTM is zero-lag cross-correlation of the source wavefield and the receiver wavefield (Etgen, 1986). The images are given by +, = ! ,, , - ., , − -
7
(3)
152 153
where is time, , and . denote the forward-propagating source wavefield and backward-propagating receiver wavefield at spatial coordinate , , and +,
154
represents the imaging result at , . The final migration image is obtained by stacking
155
the images for each event.
156
The major drawback of conventional RTM is the low-frequency noise caused by
157
cross-correlation of unwanted wavefield, which adversely affects the imaging resolution.
158
Thus, the suppression of low-frequency noise is an active research area in RTM and
159
many methods have been proposed (Baysal et al., 1984; Loewenthal et al., 1987;
160
Guittonet al., 2006; Zhang and Sun, 2008; Yoon and Marfurt, 2006). In this study, we
161
adopted the wavefield-separation method (Liu and Zhang, 2007, 2011; Xiao and Leaney,
162
2010), which applies the cross-correlation imaging condition to source wavefield and
163
receiver wavefield at different directions, as follows,
164
Ix, z = ! , ± , , - .∓ , , − - -
165
Here, the superscripts “+” and “-” denote the different directions along the horizontal or
166
vertical directions.
(4)
167 168 169
3. Synthetic test of downhole microseismic migration method For simplicity, we assume that microseismic sources, downhole sensors and
170
fractures are located in a 2D plane and accordingly design a 2D acoustic model for
171
testing the microseismic migration method introduced in the above section. Furthermore,
8
172
we assume that both source locations and velocity model are known in the numerical
173
testing of the method.
174
Fig. 1 shows the configuration of microseismic events, downhole sensors and an
175
oval-shaped fracture. The fracture is located at X=1000 m and is assumed to be induced
176
by a previous hydraulic fracturing stage. Thirty-six sensors are installed in the monitoring
177
well at X=445 m. The sensors range in depth from Z=975 m to Z=2025 m with a spacing
178
of 30 m. There are 8 microseismic events, which are assumed to be induced by a later
179
fracking stage between the fracture and the monitoring well. These microseismic events
180
are randomly distributed in the fracturing zone.
181
We used a 40 Hz Ricker wavelet as the source to model synthetic waveforms based
182
on the seismic acoustic equation. The common-event recordings are 2 minutes long. A
183
10th-order accuracy central difference stencil is applied to solve the acoustic equation
184
with the perfectly matched layer (PML) boundary condition. From the snapshot of
185
forward-propagating source wavefield at t = 0.1786 s (Fig. 2a) and the recordings (Fig.
186
2b), the scattered waves from the fractured zone can be clearly seen. For microseismic
187
migration, the background velocity model is used. The respective RTM imaging results
188
for 8 individual microseismic events are shown in Fig. 3. All images clearly reveal that
189
there exists an impedance anomaly around X=1000 m where the oval-shaped fracture is
190
located. However, the images are different for different events, indicating that the
191
contribution to imaging is different for each event because of their different locations
192
relative to the fractured zone. By stacking the normalized RTM imaging result from each 9
193
individual event, the fractured zone is better imaged due to the better data coverage (Fig.
194
4). The image of the fracture shape becomes clearer and the upper tip of the fracture is
195
clearly imaged. In comparison, the lower part of the fracture is not imaged as well as the
196
upper part. This is because there are observations from only 8 events distributed in a
197
limited aperture and the sensors are located relatively higher than the fracture. As a result,
198
sensors in the monitoring well cannot receive enough wavefield that contains information
199
related to the lower part as well as the bottom tip of the fracture.
200
There are clear artifacts in the RTM imaging result for each event due to the limited
201
spatial coverage of wavefield by using only one event (Fig. 3), and the artifacts can be
202
suppressed by using 8 events (Fig. 4). There are also some artifacts appearing to the west
203
of the monitoring well (Fig.3 and Fig. 4). This kind of artifacts has also been noted when
204
applying the RTM method to vertical seismic profile system, which is mostly caused by
205
the limited spatial coverage of the observation system and can be alleviated by the
206
increase of sensor spatial coverage (Sun et al., 2011).
207
To further test how the sensor spatial coverage affects the downhole microseismic
208
RTM imaging result, we shift the sensors both downward and upward. For the sensors
209
located from Z=1575 m to Z=2625 m, the summed RTM image becomes more accurate
210
and the shape of the fracture zone becomes clearer because more wavefields coming from
211
the lower part of fracture can be received (Fig. 5a). If the sensors are shifted upward to be
212
located in the range of Z=775 m to Z=1825 m, the fracture is more poorly resolved
213
especially for the lower part of the fracture (Fig. 5b). This indicates that with a larger 10
214
aperture of sensors, the fracture can be better characterized because wavefield from
215
different directions can be received.
216
We also test how the spacing interval of borehole sensors affect the downhole
217
microseismic RTM imaging result. For comparison, we used 12 and 24 borehole sensors
218
from Z=975 m to 2025 m with an interval of 90 m and 45 m, respectively. The fracture is
219
clearly imaged from using one event or all 8 events with 12, 24 or 36 sensors (Fig. 6). It
220
is clear, however, that the fracture is slightly better imaged with more sensors. For the
221
case of using only 12 sensors, arc-shaped spatial aliasing can adversely affect the imaging
222
resolution due to the low spatial sampling of seismic wavefield. This problem can be
223
mitigated by seismic data reconstruction, such as the compressive sensing method (Tang,
224
2010).
225
In reality, the background model is not homogeneous and for the case of shale gas
226
development the model can be approximated by a layer model and the hydraulic
227
fracturing is conducted in the target shale layer. To assess how the subsurface structures
228
affect the fracture imaging, we test two models: one by adding several layers to the
229
homogeneous background (Fig. 7a), and the other by adding a normal fault in addition to
230
layers to the homogeneous background (Fig. 7c). Except for the background velocity
231
model, all the other model setup parameters are the same as those in Fig. 1. For
232
microseismic migration in this case, we again assume the background velocity model is
233
known. The tests using two models show that even with more complex structures in the
234
background model, the fracture zone can still be clearly imaged (Figure 7b, d). At the 11
235
same time the layer interfaces below the microseismic events are also imaged, especially
236
for the first layer interface (Fig. 7b, d). Furthermore, the steep normal fault to the east of
237
the fracture zone is also clearly imaged (Fig. 7d). This indicates that for the real case of
238
hydraulic fracturing, it is possible to image both fracture zone and the shale layers by
239
using microseismic events.
240
In practice, microseismic signals may be contaminated by noise. To test how noise
241
affects the imaging result, we add Gaussian distributed noise to the recordings with
242
signal-to-noise ratio (SNR) of 10 (Fig. 8a) and 5 (Fig. 8c), respectively. The normalized
243
and stacked microseismic RTM images with all 8 microseismic events show that even
244
with noise the fracture zone can still be clearly imaged although migration images are
245
contaminated to some degree (Fig. 8b, d). This demonstrates that the
246
wavefield-separation imaging condition method used in microseismic RTM behaves well
247
in suppressing the imaging noise.
248
It is inevitable that the microseismic event is generally associated with some
249
uncertainty in its location and origin time. However, with more advanced microseismic
250
location methods such as the double-difference location method (Waldhauser and
251
Ellsworth, 2000), the microseismic event locations are very accurate and the associated
252
uncertainties could be on the order of 10 m (Castellanos and Van, 2013; Li et al., 2013).
253
Here, we further test how location errors affect the RTM imaging results by perturbing
254
event locations randomly within 10 meters (Fig. 9a), 20 meters (Fig. 9c) and 50 meters
255
(Fig. 9e), respectively. Fig. 9b, Fig. 9d and Fig. 9f show the normalized and stacked RTM 12
256
images corresponding to perturbed event locations in Fig. 9a, Fig. 9c and Fig. 9e,
257
respectively. Compared with the image using true event locations shown in Fig. 4, the
258
RTM image is almost not affected by randomly perturbing locations within 10 meters
259
(Fig. 9b). In comparison, greater location errors may adversely affect the RTM result and
260
it can be seen that the fracture structures are slightly distorted by randomly perturbing
261
locations within 20 meters (Fig. 9d). If we further perturb event locations randomly
262
within 50 m on the order of the half wavelength (Fig. 9e), the fracture zone is more
263
distorted but it is still clearly recognizable in the migrated image (Fig. 9f).
264 265 266
4. Discussion and conclusions In this study, we propose to take advantage of microseismic scattered waves and
267
apply the RTM method to directly image the fracture zone by hydraulic fracturing.
268
Unlike the conventional location-based fracture characterization method that uses
269
microseismic locations as a proxy for determining the fracture distribution, the
270
microseismic migration method that we propose to use in this study does not need to
271
detect and use all microseismic events to accurately characterize the fracture zone. As
272
long as a few microseismic events with high SNR are detected and located, they can be
273
used to characterize the fractured zone. The numerical tests showed that with only 8
274
microseismic events, we are able to clearly characterize the fracture zone.
275 276
For microseismic monitoring, it is important to estimate the SRV. More importantly, it is better to estimate the effective SRV where the proppants exist and fractures are 13
277
conductively connected (Mayerhofer et al., 2010; Johri et al., 2013). However, for the
278
zone where the microseismic events are detected and located, it does not necessarily
279
indicate the existence of proppants and opened fractures there. Therefore, it has been a
280
challenge to use microseismic monitoring to characterize where proppants migrate. With
281
the microseismic migration method, it is possible to differentiate the region with or
282
without proppants and opened fractures. Even if the region is associated with
283
microseismic events, if it does not show strong scattering effect the fractures may be
284
closed due to the lack of support by proppants.
285
We showed that for different sensor positions relative to the fracture zone, the RTM
286
imaging result is different. The monitoring system having better spatial sensor coverage
287
can better image the fracture zone. Therefore, if it is possible, the sensors should be
288
installed around the same depth zone of fracking. It is also noted that different sensor
289
spacing intervals also affect the imaging result and for the case of large spacing interval
290
the spatial aliasing caused by low spatial sampling of the wavefield will impart the
291
imaging resolution. In addition to imaging the fracture zone, from different synthetic tests
292
it is also noted that the backward-propagating wavefield is focused at the source locations
293
(Fig. 3).
294
In real cases, the microseismic signals are contaminated with noise to some degrees.
295
The synthetic tests by adding different levels of noise to microseismic recordings show
296
that the wavefield-separation imaging condition method used in microseismic RTM can
297
satisfactorily deal with the noisy data to recover the fractured zone. We also test how 14
298
microseismic event location errors affect the RTM imaging results and it shows that the
299
fracture zone can be biased to some degree but can still be clearly imaged even when the
300
event locations are randomly perturbed within 50 meters, approximately on the order of
301
half wavelength. With more advanced microseismic location methods available, the
302
location errors could well be on the order of 10s of meters. The uncertainty in the origin
303
time bears a similar effect on the RTM imaging result as location errors. It is also
304
possible that different microseismic events may have different source time functions.
305
However, due to the nature of high frequency and limited bandwidth for microseismic
306
signals, the corresponding source time functions can be well approximated by a smooth
307
ramp function with a very short duration (Li et al., 2011a, b). In addition, most of
308
microseismic events from the same fracking stage generally have similar focal
309
mechanisms. Therefore, it is expected that source signatures of microseismic events
310
would not greatly affect the RTM image.
311
In these tests, we assume the background velocity model used for microseismic
312
migration is known. In real cases, the velocity model can be well estimated from sonic
313
well logging and calibrated by perforation shots (Zhang et al., 2016). Furthermore, it can
314
also be estimated by microseismic velocity tomography using arrival times from
315
microseismic events (Zhang et al., 2009). In this study, we only use the reverse time
316
migration method to the scattered P waves, but the method can also be adapted to use
317
scattered S waves based on the elastic wave equation. We emphasize that the goal of this
318
study is to show the concept of microseismic RTM to image the fractures. For the 15
319
application of this method to the real microseismic data, we need to adapt the
320
microseismic RTM to the 3D case. One way to do this is to use the 3D acoustic equation
321
and modify the equations (1) and (2) accordingly to realize the 3D microseismic RTM.
322
Another way is to apply the 2D RTM to each microseismic event that forms a 2D plane
323
with the vertical monitoring array. By stacking RTM images from different microseismic
324
events, we may be able to obtain the 3D RTM image. To better take care of uncertainties
325
in the origin time, location and source signature of microseismic event, we could apply
326
the elastic reverse time migration imaging with both microseismic scattered P and S
327
waves (Xiao and Leaney, 2010). When combined with microseismic locations and focal
328
mechanisms, microseismic migration can help us to better understand the fracture
329
network, and thus to better optimize the fracturing design and to enhance the oil/gas
330
recovery.
331 332 333
Acknowledgements This research is supported by Natural Science Foundation of China under the grant
334
No. 41274055. We are grateful for constructive comments from two anonymous
335
reviewers that are helpful for improving this paper.
336
References
337
Baysal E, Kosloff D D, Sherwood J W C. Reverse time migration[J]. Geophysics, 1983,
338
48(11): 1514-1524. http://dx.doi.org/10.1190/1.1441434
16
339 340 341 342 343
Baysal E, Kosloff D D, Sherwood J W C. A two-way non reflecting wave equation: Geophysics, 1984, 49, 132-141. http://dx.doi.org/10.1190/1.1441644 Castellanos F, Van der Baan M. Microseismic event locations using the double-difference algorithm[J]. CSEG Recorder, 2013, 38(3): 26-37. Chambers K, Kendall J, Brandsberg-Dahl S, et al. Testing the ability of surface arrays to
344
monitor microseismic activity[J]. Geophysical Prospecting, 2010, 58(5): 821-830.
345
http://dx.doi.org/10.1111/j.1365-2478.2010.00893.x
346
Chavarria J A, Malin P, Catchings R D, et al. A look inside the San Andreas fault at
347
Parkfield through vertical seismic profiling[J]. Science, 2003, 302(5651): 1746-1748.
348
http://dx.doi.org/10.1126/science.1090711
349
Daku B L F, Salt J E, Sha L. An algorithm for locating microseismic events[C]. Electrical
350
and Computer Engineering, 2004. Canadian Conference on. IEEE, 2004, 4:
351
2311-2314. http://dx.doi.org/10.1109/CCECE.2004.1347708
352
Das I, Zoback M D. Long-period, long-duration seismic events during hydraulic fracture
353
stimulation of a shale gas reservoir[J]. SEG Technical Program Expanded Abstracts
354
2011: pp. 1473-1477. http://dx.doi.org/10.1190/1.3627480
355
Drew J, Leslie H, Armstrong P, et al. Automated Microseismic Event Detection and
356
Location by Continuous Spatial Mapping[C]. SPE Annual Technical Conference and
357
Exhibition. 2005. http://dx.doi.org/10.2118/95513-MS
358
Eisner L, Hulsey B J, Duncan P, et al. Comparison of surface and borehole locations of
359
induced seismicity[J]. Geophysical Prospecting, 2010, 58(5): 809-820.
360
http://dx.doi.org/10.1111/j.1365-2478.2010.00867.x
361
Gharti H N, Oye V, Kühn D, et al. Simultaneous microearthquake location and
362
moment-tensor estimation using time-reversal imaging[C]. 2011 SEG Annual
363
Meeting. Society of Exploration Geophysicists, 2011.
364
http://dx.doi.org/10.1190/1.3627516
17
365
Guitton A, Kaelin B, Biondi B. Least-square attenuation of reverse time migration
366
artifacts[C]. SEG Technical Program Expanded Abstracts 2006, pp. 2348-2352.
367
http://dx.doi.org/10.1190/1.2370005
368 369 370
Hill N R. Gaussian beam migration[J]. Geophysics, 1990, 55(11): 1416-1428. http://dx.doi.org/10.1190/1.1442788 Johri M, Zoback M D. The Evolution of Stimulated Reservoir Volume during Hydraulic
371
Stimulation of Shale Gas Formations. Unconventional Resources Technology
372
Conference, Denver, Colorado, 12-14 August 2013: pp. 1661-1671.
373
http://dx.doi.org/10.1190/urtec2013-170
374
Kao H, Shan S J. The source-scanning algorithm: Mapping the distribution of seismic
375
sources in time and space[J]. Geophysical Journal International, 2004, 157(2):
376
589-594. http://dx.doi.org/10.1111/j.1365-246X.2004.02276.x
377
Li J, Sadi Kuleli H, Zhang H, et al. Focal mechanism determination of induced
378
microearthquakes in an oil field using full waveforms from shallow and deep seismic
379
networks[J]. Geophysics, 2011a, 76(6): WC87-WC101. http://dx.doi.org/doi:
380
10.1190/geo2011-0030.1
381
Li J, Zhang H, Kuleli H S, et al. Focal mechanism determination using high-frequency
382
waveform matching and its application to small magnitude induced earthquakes[J].
383
Geophysical Journal International, 2011b, 184(3): 1261-1274.
384
http://dx.doi.org/10.1111/j.1365-246X.2010.04903.x
385
Li J, Zhang H, Rodi W L, et al. Joint microseismic location and anisotropic tomography
386
using differential arrival times and differential backazimuths[J]. Geophysical Journal
387
International, 2013, 195(3): 1917-1931. http://dx.doi.org/10.1093/gji/ggt358
388
Liu F, Zhang G, Morton S, et al. Reverse-time migration using one-way wavefield
389
imaging condition, 77th Annual International Meeting, SEG, 2007, pp. 2170-2174.
390
http://dx.doi.org/10.1190/1.2792917
18
391
Liu, F, Zhang G, Morton S, et al. An effective imaging condition for reverse-time
392
migration using wavefield decomposition: Geophysics, 2011, 76, s29-s39.
393
http://dx.doi.org/10.1190/1.3533914
394
Loewenthal D, Stoffa P A, Faria E L. Suppressing the unwanted reflections of the full
395
wave equation: Geophysics, 1987, 52, 1007-1012.
396
http://dx.doi.org/10.1190/1.1442352
397
Maxwell S C, Rutledge J, Jones R, et al. Petroleum reservoir characterization using
398
downhole microseismic monitoring[J]. Geophysics, 2010, 75(5): 75A129-75A137.
399
http://dx.doi.org/10.1190/1.3477966
400
Maxwell S C, Urbancic T I. The role of passive microseismic monitoring in the
401
instrumented oil field[J]. The Leading Edge, 2001, 20(6): 636-639.
402
http://dx.doi.org/10.1190/1.1439012
403 404 405
Maxwell S. Unintentional seismicity induced by hydraulic fracturing[J]. Canadian Society of Exploration Geophysics Recorder, 2013, 38(8): 40-49. Mayerhofer M, Lolon E, Warpinski N, et al. What Is Stimulated Reservoir Volume?[J].
406
SPE Production & Operations, 2010, 25(1): 89-98.
407
http://dx.doi.org/10.2118/119890-PA
408
Reshetnikov A, Buske S, Shapiro S A. Seismic imaging using microseismic events:
409
Results from the San Andreas Fault System at SAFOD[J]. Journal of Geophysical
410
Research: Solid Earth (1978–2012), 2010a, 115(B12).
411
http://dx.doi.org/10.1029/2009JB007049
412
Reshetnikov A, Kummerow J, Buske S, et al. Microseismic imaging from a single
413
geophone: KTB[C]//2010 SEG Annual Meeting. Society of Exploration
414
Geophysicists, 2010b. http://dx.doi.org/10.1190/1.3513252
415 416
Schneider W A. Integral formulation for migration in two and three dimensions[J]. Geophysics, 1978, 43(1): 49-76. http://dx.doi.org/10.1190/1.1440828
19
417
Sicking C, Vermilye J, Geiser P, et al. Permeability field imaging from microseismic[J].
418
SEG Technical Program Expanded Abstracts 2012: pp. 1-5.
419
http://dx.doi.org/10.1190/segam2012-1383.1
420
Sun W B, Sun Z D, Zhu X H. Amplitude preserved VSP reverse time migration for
421
angle-domain CIGs extraction[J]. Applied Geophysics, 2011, 8(2): 141-149.
422
http://dx.doi.org/10.1007/s11770-011-0277-1
423 424 425
Tang G. Seismic data reconstruction and denoising based on compressive sensing and sparse representation[J]. Beijing: Tsinghua University, 2010 Urbancic T I J, Shumila V J, Rutledge J T J, et al. Determining hydraulic fracture
426
behavior using microseismicity[C]. Vail Rocks 1999, The 37th US Symposium on
427
Rock Mechanics (USRMS). American Rock Mechanics Association, 1999.
428
Waldhauser F, Ellsworth W L. A double-difference earthquake location algorithm:
429
Method and application to the northern Hayward fault, California[J]. Bulletin of the
430
Seismological Society of America, 2000, 90(6): 1353-1368.
431
http://dx.doi.org/10.1785/0120000006
432
Xiao X, Leaney W S. Local vertical seismic profiling (VSP) elastic reverse-time
433
migration and migration resolution: Salt-flank imaging with transmitted P-to-S
434
waves[J]. Geophysics, 2010, 75(2): S35-S49. http://dx.doi.org/10.1190/1.3309460
435
Yang Y, Zoback M D. The role of preexisting fractures and faults during multistage
436
hydraulic fracturing in the Bakken Formation[J]. Interpretation, 2014, 2(3):
437
SG25-SG39. http://dx.doi.org/10.1190/INT-2013-0158.1
438 439 440
Yoon K, Marfurt K J. Reverse-time migration using the poynting vector [J]: Exploration Geophysics, 2006, 59, 102-107. http://dx.doi.org/10.1071/EG06102 Zhang H, Wang P, van der Hilst R D, et al. Three-dimensional passive seismic waveform
441
imaging around the SAFOD site, California, using the generalized Radon
442
transform[J]. Geophysical Research Letters, 2009, 36(23).
443
http://dx.doi.org/10.1029/2009gl040372 20
444
Zhang, H., S. Sarkar, M.N. Toksoz, H.S. Kuleli, and F. Al-Kindy, Passive seismic
445
tomography using induced seismicity at a petroleum field in Oman, Geophysics, 2009,
446
74(6), WCB67, doi:10.1190/1.3253059.
447
Zhang, Jiewen, H. Zhang, Y. Zhang, and Q. Liu, Fast one-dimensional velocity model
448
determination for microseismic monitoring using station-pair differential arrival times
449
based on the differential evolution method, Physics of the Earth and Planetary
450
Interiors, in revision 2016.
451
Zhang Y, Sun J. Practical issues of reverse time migration-true-amplitude gathers, noise
452
removal and harmonic-source encoding[C]. 70th EAGE Conference & Exhibition.
453
2008. http://dx.doi.org/10.3997/2214-4609.20147708
454
Zhebel O, Eisner L. Simultaneous microseismic event localization and source mechanism
455
determination[J]. Geophysics, 2014, 80(1): KS1-KS9.
456
http://dx.doi.org/10.1190/segam2012-1033.1
457
Zoback M D, Kohli A, Das I, et al. The importance of slow slip on faults during
458
hydraulic fracturing stimulation of shale gas reservoirs[C]. SPE Americas
459
Unconventional Resources Conference. Society of Petroleum Engineers, 2012.
460
http://dx.doi.org/10.2118/155476-MS
461 462 463 464 465 466 467 468 469 470 21
471
Figure Captions
472 473
Fig. 1. 2D P-wave velocity model with a fracture zone at X=1000 m. The background
474
velocity is 4200 m/s and the fracture zone is associated with a low velocity anomaly of 15%
475
lower than the background velocity. Monitoring well is located at X=445 m with 36
476
sensors (denoted by triangles) ranging in depth from Z=975 m to Z=2025 m with a
477
spacing of 30 m. Microseismic events are (denoted by asterisks) located at (850 m, 2120
478
m), (850 m, 1975 m), (850 m, 2260 m), (862 m, 2048 m), (862 m, 2170 m), (830 m, 2107
479
m), (843 m, 2041 m), and (842 m, 2235 m), respectively.
480 481
Fig. 2. (a) Snapshot of source wavefield at t=0.1786 s. (b) Common-shot recordings for
482
the event located at (850 m, 1975 m) and the sample interval is 0.595 ms. The model
483
setup is shown in Fig. 1.
484 485
Fig. 3. RTM imaging result from 8 individual microseismic events (denoted by asterisks)
486
for the model setup in Fig. 1.
487 488
Fig. 4. Normalized and stacked RTM imaging result from 8 microseismic events shown
489
in Fig. 3.
490 491
Fig. 5. Comparison of RTM imaging results for sensors located at different depths for the
492
model setup shown in Fig. 1. (a) From Z=1575 m to 2625 m. (b) From Z=775 m to 1825
493
m.
494 495
Fig. 6. Comparison of imaging results using (left) the individual event located at (850 m,
496
2120 m) and (right) all 8 events together in different cases of having (top) 12 sensors,
497
(middle) 24 sensors and (bottom) 36 sensors. The span of these sensors is the same in 22
498
depth but with different spacing of 90, 45 and 30 m from top to bottom. Other symbols
499
are the same as those in Fig. 3.
500 501
Fig. 7 Comparison of the RTM images when the background velocity models are
502
complex. (a) The background velocity model contains several layers with different
503
velocities compared to that shown in Fig. 1; (b) Normalized and stacked RTM image
504
from 8 microseismic events for the model setup in (a); (c) The same as (a) but with a
505
normal fault to the east of the fracture zone; (d) Normalized and stacked RTM image for
506
the model setup in (c). The model layers are denoted by dashed lines in (b) and (d).
507 508
Fig. 8 The RTM images from microseismic recordings with different levels of noise. (a)
509
The recordings from a microseismic event by adding Gaussian distributed noise with
510
SNR of 10. (b) Normalized and stacked RTM image for 8 events for the noise level in (a).
511
(c) The recordings from a microseismic event by adding Gaussian distributed noise with
512
SNR of 5. (d) Normalized and stacked RTM image for 8 events for the noise level in (c).
513
The model setup is the same as Fig. 1.
514 515
Fig. 9 Normalized and stacked RTM images (right) from 8 microseismic events when
516
their locations are randomly perturbed at different degrees (left). (a) and (b) show
517
locations randomly perturbed within 10 m and the corresponding RTM image. (c) and (d)
518
show locations randomly perturbed within 20 m and the corresponding RTM image. (e)
519
and (f) show locations randomly perturbed within 50 m and the corresponding RTM
520
image. White asterisks denote true seismic locations and red asterisks denote perturbed
521
locations.
522
23
523 524
24
525 526
25
527 528
26
529 530
27
531 532
28
533 534
29
535 536
30
537 538
31
539
32