Current developments on micro-seismic data processing

Accepted Manuscript Current developments on micro-seismic data processing Hua Wang, Meng Li, Xuefeng Shang PII:

S1875-5100(16)30103-2

DOI:

10.1016/j.jngse.2016.02.058

Reference:

JNGSE 1308

To appear in:

Journal of Natural Gas Science and Engineering

Received Date: 30 October 2015 Revised Date:

21 January 2016

Accepted Date: 25 February 2016

Please cite this article as: Wang, H., Li, M., Shang, X., Current developments on micro-seismic data processing, Journal of Natural Gas Science & Engineering (2016), doi: 10.1016/j.jngse.2016.02.058. 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.

ACCEPTED MANUSCRIPT

Current developments on micro-seismic data processing Hua Wang, Earth Resources Lab. MIT

Xuefeng Shang, Earth Resources Lab. MIT

Abstract

RI PT

Meng Li, China University of Petroleum;

SC

In the past decades, micro-seismicity, mainly induced by fluid injection (or depletion),

M AN U

has attracted considerable interest due to its ability to monitor the development of unconventional oil/gas field, such as volumetric extension of hydro-fractures in tight and shale gas fields as well as geo-hazards due to reservoir subsidence. With the deployment of portable (e.g. surface and downhole geophones) and permanent (such

TE D

as optical fiber) monitor systems, invaluable information of fluid flow path, local stress changes, geomechanical deformation can be inferred from micro-seismic data, though it is still very challenging to interpret the results without any ambiguity. In this

EP

paper, we briefly review some key ingredients of micro-seismic data processing, such as data acquisition, characteristics of micro-seismic events, event onset time

AC C

auto-picking, event location and relocation as well as source mechanism analysis. The uncertainties and issues of current processing workflow are addressed, which is critical to integrate the micro-seismicity into reservoir monitor tool. Finally potential improvement and application are discussed, such as new distributed acoustic system.

ACCEPTED MANUSCRIPT

Introduction In the past decades, fluid injection (or depletion) treatments are used to improve

RI PT

underground resource recovery, such as hydraulic fracturing in tight sand and shale oil/gas reservoirs with very low permeability (Economides and Nolte, 2000), steam

injection for the heavy-oil fields, and enhanced geothermal systems in the geothermal

SC

resource development (e.g. Miyazawa et al., 2008; Maxwell et al., 2010). Massive fluid injection (or depletion) considerably changes the subsurface stress in the

M AN U

reservoir, which results in new fracture development or pre-existing fracture reactivation. The first case study of injection induced seismicity is in a waste water injection well in Denver, USA (see Healy et al., 1968). Since then induced seismicity monitoring is increasingly popular during the production, especially as the surge of

TE D

unconventional oil/gas. For example, Lei et al. (2013) studied the fluid injection-induced seismicity in a natural gas reservoir. Martínez-Garzón et al. (2014) observed the spatial-temporal distribution pattern of induced seismicity in a

EP

geothermal reservoir. Skoumal et al. (2015a and 2015b) show case studies on induced

AC C

earthquakes induced by hydraulic fracturing and micro-seismicity induced by deep waste water injection.

The occurrence of such fluid-injection-induced seismicity strongly correlates with rock brittleness characterized by high Young’s modulus and low Poisson’s ratio (e.g. Langenbruch and Shapiro, 2015). It is also highly correlated with the treatment pressure and injection rate (e.g. Sleefe, 1995). Different rock failure criteria may

ACCEPTED MANUSCRIPT

apply. In hydraulic fracturing, the fracture tip opens due to extensile failure, and microseisms are induced close to the tip zone. The location of microseisms distribute

RI PT

along the new fracture as the fracture tip propagates, from which the dimension and geometry of new factures can be inferred (Pearson, 1981; Warpinski et al., 2001). For

pre-existing fracture/fault reactivation, Coulomb shear failure criterion usually applies,

SC

in which the orientation of the fracture roughly aligns with the maximum shear stress direction.

M AN U

In conventional production field, the micro-seismicity behaves different from above. The magnitude ranges from -3 to 4, whereas magnitude resulting from hydraulic fracturing for fracture opening is below 0. The locations are not concentrated within the reservoir layer. People speculate it is related with the reservoir properties,

Dye et al., 1999).

TE D

lithology, rock properties of the vicinity layers and pressure drop of the reservoir (e.g.

EP

Different from conventional exploration seismic (surface survey), the information of source such as the location, magnitude, fire time and focal mechanism is unknown in

AC C

micro-seismicity monitoring, which makes it more complex to process and interpret. The first and most accurate parameter inferred from micro-seismicity data is the event location. For instance, the spatial distribution of the events in hydro-fracturing roughly delineates the fracture planes and networks around the well, on which the reservoir fluid model is based (Baig and Urbancic, 2010). The event location does not always, however, coincide with the new rock fracture locations, and sometimes it is

ACCEPTED MANUSCRIPT

an indicator of pre-existing fracture/fault reactivation. Most of small induced microseismic events in the pre-existing fault/fracture slip can be approximated by a

RI PT

double-couple point source (Rutledge and Phillips, 2003). However, the hydraulic fracture opening contains significant volumetric change (ISO) or compensated linear vector dipole (CLVD) component ( ílený, 2009). Combining the location and focal

SC

mechanism, we can infer the plane orientation of factures, connectivity of fracture

network, fracture opening mechanism and the local stress field distribution around the

M AN U

facture (Adam and Urbancic, 2010). We can also estimate the optimal interval for perforation and hydraulic fracturing wells and stimulated reservoir volume (SRV), which is critical to increase the efficiency and reduce the cost of the operation (Maxwell, 2011). The application in the conventional reservoir are summarized by

TE D

Dye et al. (1999), such as identifying fault structures that can result in reservoir compartmentalization, fracture-dominated

imaging

reservoirs,

the

flow

providing

anisotropy real-time

4-D

associated monitoring

with of

EP

fluid-pressure-front movement, assisting in targeting new injector/producer wells, identifying potential wellbore instability and so on.

AC C

Many typical case studies ve already conducted in recent years. Stabile et al. (2014) reveals a dipping fault using induced seismicity. Block et al. (2015) use the precise relative location results of induced seismicity to constraint the subsurface geological structure. Improte et al. (2015) consider the induced seismicity rate correlates with fluid injection curves and temporal variations of elastic and anisotropic parameters. Wei et al. (2015) studies the 2012 Brawley swarm triggered by injection-induced

ACCEPTED MANUSCRIPT

aseismic slip. Will et al. (2014) discover the relationship between micro-seismic event occurrence and subsurface geology in multiple domains and at multiple scales.

RI PT

Hornbach et al. (2015) summarize the causal factors for seismicity along a mapped ancient fault system near Azle, Texas. Rutledge et al. (2015) studied the

source-mechanism trend for the interpretation of the microseismic shearing driven by

SC

hydraulic fracturing.

M AN U

In this paper, we briefly review recent progress on the data acquisition, processing and applications in micro-seismic monitoring system and workflow.

Micro-seismic Data Acquisition

TE D

A micro-seismic monitoring system consists of sensors, data acquisition, data storage and transmission, data processing and imaging, as well as real time feedback. The following sections will give the details of those different parts.

EP

Monitoring Sensors and Data acquisition

AC C

Various sensors for micro-seismic monitoring are available, including surface sensors (portable earthquake station with short period, geophone with high frequency and low-sensitivity (Duancan and Eisner, 2010)) and downhole sensors, similar to sensors in vertical seismic profiling (VSP) but in higher frequency band (~ 2 kHz). The surface sensors are placed with dense spacing and good coverage with a low cost. But the signal-to-noise ratio (SNR) is relatively low due to large distance from the

ACCEPTED MANUSCRIPT

underground events and near surface noise contamination. Higher frequency component in the signal are attenuated by the formation as well. Nevertheless, the

RI PT

surface sensors work well in an exposed basement area where the near surface effect is small. Figure 1 shows an example of a surface monitoring system in the

southwestern China where the base rock is less than half meter depth beneath the

SC

surface (Zeng, et al., 2014). The sensor is specified for local earthquake observations, consisting of three component measurements of ground velocity, GPS (Global

M AN U

Position System) antenna, battery and data collection system. The origin point in Figrue 1 (right) shows the location of a shale gas hydraulic fracturing wellhead, which is surrounded by the surface monitoring system. On the other hand, the downhole sensor is usually placed in a nearby well (about 300 meters far away from the

TE D

treatment well) and can record weak events with high SNR. However, the monitoring well is expensive and is not accessible in a new field. The data coverage, especially the azimuthal coverage, is limited, which introduces uncertainties to the event

EP

location and source mechanism. Figure 2 shows the schematic of a down-hole 3C sonde and an example of the downhole monitoring system. In the example shown in

AC C

Figure 2(right), there are 15 hydraulic fracturing stages in the horizontal treatment well, in which the plug and perforation for different stages are labeled by different colors. A monitoring well with 12 downhole sensors (30 m interval) is around 300 meters away from the treatment well.

RI PT

ACCEPTED MANUSCRIPT

SC

Figure 1 Surface monitoring system of micro-seismicity: (Left) A surface monitoring sensor system with a short-period 3C sensor, GPS antenna, battery and data collection system;

(Right)

an example for the placement of the surface monitoring system around a shale gas hydraulic

M AN U

fracturing well in southwest China. Horizontal and vertical axes are the east-west and north-south direction. The original point of the map is the well site. The triangles are the locations of surface sensors and the line in the middle of figure is the horizontal well trajectory.

EP

TE D

(Figure 1b from Zeng et al., 2014)

AC C

Figure 2 Schematic of a down-hole 3C sonde and (left) and an example for the placement (right)

Data storage and transmission

Ideally, micro-seismic monitoring system transmits the data to the operator (or processing center) in real time. In the hydraulic fracturing operation, the operator benefits from the real time micro-seismic location in treatment quality control.

ACCEPTED MANUSCRIPT

However, the limitation of the field condition makes the real time transmission challenging in some situations, and the data is stored in disk and analyzed afterwards.

RI PT

Characteristics of the micro-seismic signal The signals of induced microseismicity are dominated by shear motion with larger

amplitude of shear (S) wave than that of compression (P) wave. The frequency

SC

usually ranges from a few Hz to thousands Hz depending on the distance between the

M AN U

sensors and the treatment well. We will show some typical events during the micro-seismic monitoring to understand the features of the micro-seismic signal.

Figure 3 shows a downhole monitoring example of the perforation shots which are

TE D

used to calibrate the orientations of 3-C geophones during micro-seismic monitoring (from Figure 3 in Eaton et al., 2013). The left three figures are the frequency spectrum (top), trace (middle), and time-frequency displays (bottom) for the event in a sensor

EP

near the perforation shot. The right figures are from a sensor far from the perforation shot. It is clear that the signals from perforation shots are dominated by P arrivals

AC C

because the perforation treatments generate volumetric explosive sources. As the distance from the shot increases, the noise become larger and the high frequency component of the signal is attenuated, which makes the shot signal barely detected by surface monitoring sensors.

RI PT

ACCEPTED MANUSCRIPT

SC

Figure 3 Examples of perforation shots, which are generally dominated by P arrivals. Upper

panels: signal and noise spectra. Lower panels: recorded waveforms and their corresponding

M AN U

spectrograms obtained by using a short-time Fourier transform (From Figure 3 in Eaton et al., 2013)

Figure 4 is an example for a micro-seismic event (Mw = 0.8) induced by fluid

TE D

injection in a geothermal field. The magnitude of the event is small and buried in strong ambient noise (upper figure). The signals with coherence stand out in 32s and

EP

247s (frequency above 2 Hz) based on time-frequency analysis (middle figure). After a band-pass filter from 2 to 20 Hz, we can get a clear signal (bottom figure). This is a

AC C

good example how to extract weak events in noisy data (low SNR) by using time-frequency analysis.

SC

RI PT

ACCEPTED MANUSCRIPT

M AN U

Figure 4 A micro-seismic event (Mw = 0.8) induced by fluid injection in a geothermal field. Top: the raw data; Middle: time-frequency display with a short time Fourier transform; Bottom: the filtered data with a band-pass filter from 2 to 20 Hz.

Figure 5 shows a typical event and the spectrograms recorded by a borehole network

TE D

with depth between 750 m to 1250 m from surface (From Figure 4 in Li et al., 2011b). We can find the frequency of the recorded waveform is high where the significant

AC C

EP

energy ranges from 15 to 35 Hz.

ACCEPTED MANUSCRIPT

Figure 5 The seismograms of a typical event recorded by the borehole network. The filtered seismograms (15-35 Hz) are in the left column; the original seismograms are in the middle; the spectrograms of the original seismograms are in the right. The borehole data are dispersive, i.e., higher frequency contents arrive later as the energy is trapped within layers and propagates as

RI PT

guided waves. (Figure 4 in Li et al., 2011b)

SC

Data Processing Workflow

A standard work flow of micro-seismic processing is listed as below. Firstly, the

M AN U

orientation of three component geophone data and velocity model are calibrated by check shots (or perforation shots). After that, the noise is suppressed by filtering and time-frequency to enhance the signal-to-noise ratio. Subsequently, the event auto-detector, such as STA/LTA, AIC and Karlman filter, marks the onset time of

TE D

potential micro-seismic events. The false alarms are then removed by QA/QC process. The detected events can be located by either ray or full waveform methods based on the velocity model and onset time. Double difference (DD) method can further

EP

minimize the influence on the event location of velocity model uncertainty. Source mechanism can be derived by moment tensor inversion as well, if acquisition

AC C

coverage is good or the shear wave arrival time is clear. Finally, the micro-seismic image is interpreted as the distribution of fluid flow and induced fractures.

ACCEPTED MANUSCRIPT

Micro-seismic Data Input

Check shot Calibration

Orientation Calibration

Polarization Analysis of 3C Data of Check shots

Noise Removal

Filtering & Time-Frequency

Micro-seismic Events Detection

STA/LTA, AIC, Karlman

Micro-seismic Imaging and Interpretation

SC

M AN U

Event Location

RI PT

Velocity Model Calibration

Grid Searching, RTM & DD

TE D

Figure 6 A standard workflow for Micro-seismic processing

EP

Micro-seismic Event Detection

AC C

The micro-seismic data are continuously recorded, so the first step is the effective event detection and picking from the background noise. Various methods for detecting earthquakes in seismology have been introduced into micro-seismic processing. However, not all of them can be directly used due to the distinctive characteristics of the micro-seismic signals. The energy ratio, Akaike information criteria (AIC), neural network, fractal theory, Karlman filter (Kalman, 1960) and their combinations are currently applied to detect and pick up the events.

ACCEPTED MANUSCRIPT

The short term average/long term average (STA/LTA) method, one of the most

RI PT

classical energy ratio methods, was first addressed by Allen in 1978 (Allen,1978). By calculating the ratio of the characteristic functions in a long and short time window, we can extract the events and their corresponding arrival time from the noisy data if

SC

the ratio is larger than a pre-defined threshold. The method is stable and suitable for

real-time automatic processing of the massive data. However, the determination of the

M AN U

optimized window size and threshold is the key for the success of this method. A shorter window is more sensitive to small abrupt changes in the waveform, allowing us to pick the weak signals, but is also easier to be influenced by strong noise and vice versa. Similarly, a small threshold can detect signals in low S/N situations but the

TE D

noise might also be flagged as a true event. To minimize the influence of the parameters, many modifications are introduced on characteristic function combined with polarization analysis in multi scale (Kanwaldip and Dowla, 1997), adaptively

EP

choosing the window size based on the wavelet transform as well as determining threshold based on noise level and additional processing window (Chen and Stewart,

AC C

2006).

Equation 1 shows the characteristic function CF, where X refers to the amplitude of received waveforms and i is the data point. Then, the short and long term average and their ratio can be calculated as shown in equation 2, where Nsta and NLta is the number of data points in short and long processing window.

ACCEPTED MANUSCRIPT

CF ( i ) = X ( i ) − X ( i − 1) X ( i + 1) 2

STA (i ) =

(1)

1 Nsta ∑ CF (i) Nsta i =1

1 Nlta ∑ CF (i) Nlta i =1 STA(i ) Ratio(i ) = LTA(i ) LTA(i ) =

RI PT

(2)

Here is an example to illustrate the features of the STA/LTA method, as shown in

SC

Figure 7. Figure 7a shows the original waveform where the P wave and S wave are

clearly observed. Figure 7b shows the STA (blue) and LTA (red) and Figure 7c shows

M AN U

the STA/LTA ratio (blue) as well as the threshold (red). The threshold is usually set as the same level as the SNR of the input trace, which is estimated by the RMS of pre-signal data. Based on Figure 7, it is obvious that the peak of STA/LTA ratio is larger than the pre-defined threshold, which suggested that the input signal contains For noise-free seismograms, the maximum value of the

TE D

the effective events.

numerical derivative of the STA/LTA ratio is close to the onset.

EP

Amplitude

1000

(a)

0

−1000

0

0.5

1

1.5

2

2.5

3

Time(s)

3.5 −3

x 10

4

x 10

STA,LTA

AC C

5

0

0

STA

(b)

LTA

0.5

1

1.5

2

2.5

3

Time(s)

3.5 −3

x 10

STA/LTA

2000 1000 0

0

(c)

STA/LTA Threshold

0.5

1

1.5

2 Time(s)

2.5

3

3.5 −3

x 10

Figure 7 Micro-seismic events picking by STA/LTA: (a) Original waveform; (b) STA (blue) and LTA (red) results; (c) STA/LTA ratio (blue) and the threshold (red).

ACCEPTED MANUSCRIPT

Autoregressive-Akaike information criteria (AR-AIC) method is another effective

RI PT

method to detect P and S phases. This method assumes that the recorded waveform segments considered as an autoregressive (AR) process and the segments before and after the onset are assumed different (Sleeman and Eck, 1999). As shown in equation

SC

3, where A is the AIC value, k is the merging point for two-interval model for seismogram with length of N, M is the order of an AR process fitting the data, C2 is a constant and σ 1,max andσ 2,max indicate the variance of the seismogram in the two 2

M AN U

2

intervals. The sample point that makes the AIC value reaches its minimum is considered as the arrival time. The advantage of this method is that it provides more precise picking results compared to energy ratio method like STA/LTA when the S/N

TE D

is high. However, it is found that AIC is not reliable for the data with strong noise. In addition, AIC would always get a pick no matter whether it is signal or not. Therefore, a precise window length is needed. Zhang and Thurber (2003) assumes that the

EP

primary features such as P and S onsets are retained and other features such as scatter arrivals and noise decay quickly after wavelet decomposition. Therefore, they first use

AC C

discrete wavelet transform to decompose the original waveform into several scales (several wavelet coefficients). Then a threshold is applied to each wavelet coefficient to suppress the random noise. After that, the AIC picker is applied to each wavelet coefficient to get pickings at each scale. If the AIC pickings are consistent (the time difference of each picking is within a pre-defined limit), an effective event is claimed, otherwise, no event exist in the original seismogram.

ACCEPTED MANUSCRIPT

2 2 A(k ) = (k − M ) log(σ 1,max ) + ( N − M − k ) log(σ 2,max ) + C2

(3)

Figure 8a and Figure 9a show the noise and signal segment of a recorded waveform.

RI PT

Figure 8 (b-d) and Figure 9 (b-d) show the AIC picking results at different levels for noise and signal segment, respectively. It is obvious that the AIC picking results in different scales are not consistent in noise segment while those are very similar in

SC

signal segment. The red line in Figure 9 (a) represents the manual picking, which is

M AN U

close to multi-scale AIC picker.

Some other methods for event picking are also in development. Song et al. (2010) introduced a matched filter method to detect the events and pick the arrival time, which is an array based waveform correlation approach. Although some modifications

TE D

and revisions have been introduced in detection algorithms, precise and automatic detection and pick of effective micro-seismic events is still changeling.

EP

Since micro-seismic measurements are often mixed with strong environmental noise, picking results of direct application of these methods seem to be not reliable and Therefore,

AC C

stable.

application

of

reasonable

de-noising

algorithm

before

micro-seismic processing becomes critical especially in surface micro-seismic monitoring where S/N is commonly lower than 1 dB. In addition, most picking algorithms can only extract the first arrival of P wave. The current strategy to identify S wave is to reapply the detection method in a short time interval after P wave first arrival. However, this strategy needs further analysis such as hologram analysis to

ACCEPTED MANUSCRIPT

judge whether the second detected phase belongs to S wave or another P wave. Therefore, development of automatic detection and pick of both P wave and S wave

0.05 0 −0.05

(a)

−0.1 0

0.5

1

1.5

2

Time(s)

AIC

x 10

−4.213

AIC values at scale 1 (b)

−4.214

0

0.5

1

1.5 Time(s)

4

AIC

−4.196 −4.198

2

M AN U

x 10

2.5

AIC values at scale 2

(c) −4.2

2.5

SC

4

−4.212

RI PT

Amplitude

algorithm is urgent.

0

0.5

1

1.5

2

2.5

2

2.5

Time(s)

4

AIC

−4.482

x 10

−4.484

AIC values at scale 3

(d) −4.486

0

0.5

1

1.5

Time(s)

TE D

Figure 8 Multi-scale AIC picker of a noisy data: (a) noisy data; (b) to(d) are AIC picking results at levels of 1 to 3, respectively

0

EP

Amplitude

1 0.5

(a)

−0.5

1.5

2

2.5

3

3.5

4

Time(s)

4

AIC

−4.08

x 10

−4.09

AIC values at scale 1

AC C

(b)

−4.1 1.5

2

2.5

3

3.5

4

Time(s)

4

AIC

−3.4

x 10

−3.6

AIC values at scale 2

(c)

−3.8 1.5

2

2.5

3

3.5

3

3.5

4

Time(s)

4

AIC

−3

x 10

−3.5 AIC values at scale 3

(d) −4 1.5

2

2.5

4

Time(s)

Figure 9 Multi-scale AIC picker of a signal with good SNR: (a) data segment with good SNR; (b) to(d) are AIC picking results at levels of 1 to 3, respectively

ACCEPTED MANUSCRIPT

RI PT

Micro-seismic Event Location

The locations of micro-seismic events are of great importance and are the only way

for determining the distribution and geometry of the fractures induced by hydraulic

SC

fracturing. In the last few years, a variety of micro-seismic events location methods have been tested and developed where most of them are introduced from seismology.

M AN U

Among them, the Geiger’s method (Geiger, 1912) is the most classical linear inversion method for earthquake location. Subsequently, many modifications and revisions have been introduced to improve the performance of Geiger’s method such as the singular value decomposition and damp least square (Lienert et al., 1986)

TE D

method to increase the stability of linear equations as well as the conjugate gradient method or non-linear optimization methods like simplex (Prugger and Gendzwill, 1988) and global grid searching (Nelson and Vidale, 1990) to improve the efficiency

EP

in event location.

AC C

Generally, the micro-seismic event location methods can be divided into two categories. The first one is based on the first arrivals of P and/or S wave only while the other one is based on waveform amplitude. The methods based on arrival time employs some ray tracing techniques to minimize the difference between observed and calculated arrival times in a certain velocity model. The grid point which has the minimum of arrival time residual is considered as the event location. Probabilistic grid searching method (Lomax et al., 2001), conjugate gradient or any other global

ACCEPTED MANUSCRIPT

optimization algorithm such as Monte Carlo (Hammersley and Handscomb, 1967) and simulated annealing (Billings, 1994) methods can be employed to find the

Probabilistic Grid Searching Method

RI PT

optimized solution.

Probabilistic grid searching method is based on the picked arrival time, probability

SC

density function and ray tracing technique. In addition, the method includes the picking errors of different phases to obtain the location results and accuracy

M AN U

simultaneously. The velocity model is first divided into many grid points. Then, an objective function of location is built, as shown in equation 4. X refers to the grid point coordinates X = (x,y,z), N is the receiver number, tPobs and tSobs are the arrival times of picked P and S waves from received waveform, tPcal and tScal are calculated

TE D

arrival times of P and S waves from receiver to the grid point by using the ray tracing technique, N is the number of receivers, σS andσP is the standard deviations which

EP

summarize the assigned uncertainties on picked P and S arrivals. The grid point which

AC C

makes the objective function reach its minimum is assumed as the source location. N

 [(t obs − tPobs ) − (tScal − t Pcal )]2  1 L ( x) =  ∑ exp(− S ) 2 2 2 2 σ + σ N + σ σ S P S P   (4)

Collapsed grid searching algorithm is an effective strategy to find the optimized location, Figure 10 illustrates the principle of this method. Firstly, the coarse grids are divided based on the given velocity model to calculate the probability density function (PDF) of the location. Then, a certain grid with smallest arrival time residue is

ACCEPTED MANUSCRIPT

subdivided into fine grids. The dividing procedure is continued until the arrival time

M AN U

SC

RI PT

residue is less than a pre-defined threshold.

Figure 10 Collapsed grid searching principle: (A) ture probability density function (PDF) of

TE D

the seismic location with coarse grid; (B) to (F) different subdivisions with fine grids for location.

Figure 11 shows a 3D numerical example of micro-seismic event location by using grid searching method based on synthetic data. Fast marching method (FMM) is

EP

chosen as the ray tracing technique to calculate the travel time between receivers and

AC C

each grid point in forward modeling and inversion (Sethian and Popovici, 1999). Global grid searching strategy is applied to search the optimized solution to the objective function (equation 4). Figure 11 (a) shows a four-layered velocity model, receiver geometry (red triangles) and the locations of micro-seismic events (blue dot). The model has 200 grid points in X, Y and Z direction and the grid spacing used in FMM is 10m. Figure 11 (b) and (c) show the East - North and North-Depth view of probabilistic grid searching results whenσS =σP= 0.1s, respectively. Figure 11 (d)

ACCEPTED MANUSCRIPT

and (e) show the East - North and North-Depth view of probabilistic grid searching results whenσS =σP= 0.05s, respectively. It is obvious that fig.11 (d) and (e) has a

RI PT

smaller location radius representing a higher location accuracy than fig.11 (b) and (c). SinceσS andσP measures the accuracy of S and P wave arrival pickings, higherσS and σ P means lower location accuracy. The detailed uncertainty analysis of

200 180

North<->South(m)

160 140 120 100 80 60 40 20

1

(c)

0.9

20

0.8

40

0.7

60

0.6

80

0.5 0.4

Depth(m)

(b)

M AN U

SC

probabilistic grid searching method can be found in Leo (2009).

1

0.8

0.6

100 120

0.3

140

0.2

160

0.1

180

0.4

0.2

200

200 180

North<->South(m)

160 140 120 100 80 60

1

0.4

0.2

EP

40

100 East<->West(m)

50

100 East<->West(m)

150

200

0

1

40

0.8

60 80

0.6

100 120

0.4

140 160

0.2

180

20

50

200

20

0.8

0.6

150

(e)

Depth(m)

(d)

100 East<->West(m)

TE D

50

150

200

0

200 50

100 East<->West(m)

150

200

0

AC C

Figure 11 Micro-seismic event location using probabilistic grid searching based on FMM

and VFSA: (a) a four-layered velocity model, receiver geometry (red

traiglues) and the

locations of micro-seismic events (blue dot). Color bar represents the velocity value, which varies from 3000m/s to 4500m/s. (b) and (c) East - North and North-Depth view of probabilistic grid searching results whenσS =σP =0.1s, respectively; (d) and (e) East - North and North-Depth

view of probabilistic grid searching results whenσS =σP =0.05s, respectively. Colorbars for (b) to (e) show the normalized energy of probabilistic grid searching images from 0 to 1. The highest energy in the images corresponds to the event location.

ACCEPTED MANUSCRIPT

Intersection Method The intersection method (Zhao et al., 2008; Zhou and Zhao, 2012) is based on the

RI PT

reciprocity of wavefield. The receivers are treated as sources in this method and estimate the arrival time from these receivers to each grid point. Then we use

observed travel time difference of P wave and S wave at a certain receiver or P wave

SC

at adjacent receivers as the constraints to get the trajectory with same travel time residue. Multiple trajectories can be obtained after we apply the same processing to

location of micro-seismic event.

M AN U

each of the receivers. The intersection of all the trajectories gives us the estimated

Equation 5 shows the distribution of travel time difference of P wave and S wave. DTPScal represents the calculated travel time difference between P wave (Tp) and S

TE D

wave (TS) from the jth receiver to grid point x = (x,y,z), N is the receiver number.

DTcal PS ( x, R Pj , R Sj ) = T P ( x, R j ) − T S ( x, R j ), j = 1, 2,K N

(6)

EP

FDDT PS ( x, R jP , R Sj ) =| DTcalPS ( x, R Pj , R Sj ) − (τ P R j − τ S R j ) |

(5)

AC C

The difference between calculated travel time difference and observed travel time difference can be estimated by equation 6, where DTPScal represents the calculated

travel time difference between P wave (Tp) and S wave (TS) from the jth receiver to grid point x = (x,y,z), τ P R j and τ S R j is the picked arrival time of P and S wave, respectively. Figure 12 shows a 2D synthetic example of event location by using intersection method. Figure 12 (a) and (b) show the velocity model and traced rays in a given

ACCEPTED MANUSCRIPT

receiver geometry and a single micro-seismic event. Figure 12 (c) shows the distribution of travel time residue in the whole model of the travel time difference of P

RI PT

and S waves at the first receiver. The yellow area in the figure is the trajectory, which also represents the minimum of travel time residue between observed and calculated time difference. Figure 12 (d) shows the distribution of travel time residue and

SC

trajectory of P wave travel time difference between the first and second receiver.

Figure 12 (e) and (f) shows all the trajectories of P wave and S wave travel time

M AN U

difference at each receiver as well as P wave travel time difference at adjacent receivers. The intersections (Red Cross) in these two figures stand for the estimated location of micro-seismic events. It is found that the trajectories of P wave and S wave travel time difference are mainly constrained in Z direction while trajectories of

TE D

P wave travel time difference at adjacent receivers are more constrained in X direction. It means that the travel time differences of P and S waves has a higher accuracy in Z direction while P wave travel time difference at adjacent receivers has a higher

AC C

EP

accuracy in X direction.

ACCEPTED MANUSCRIPT

Figure 12 Micro-seismic event location using intersection method based on FMM and VFSA: (a) velocity model, red star is ture event location and red points are receivers; (b) traced rays of the model; (c) distribution of travel time residue in the model of the travel time difference of P and S waves at the top receiver. The yellow area is the trajectory of the minimum of travel time

RI PT

residue between observed and calculated time; (d) distribution of travel time residue and trajectory of P wave travel time difference between the first and second receiver; (e) and (f) all the trajectories of P wave and S wave travel time difference at each receiver as well as P wave

travel time difference at adjacent receivers. Red cross is the estimated location of micro-seismic

SC

events.

Travel-time based location methods such as grid searching and intersection method

M AN U

are easy to be implemented and stable with high efficiency. As discussed before, these methods are based on travel time only without the use of waveform amplitude, which means the inaccurate picked arrivals might greatly degrade the accuracy and reliability of location results especially in noisy data. In addition, the inaccurate

TE D

velocity model definitely introduces great errors in location results. Therefore, the locations obtained by this method can only be considered as preliminary results. It requires other relative location methods such as double difference to minimize the

EP

influence of the inaccurate velocity model for the refinement of location results.

AC C

Reverse Time Imaging Method As mentioned above, the picked arrivals are not reliable due to the strong background noise, which greatly degrades the accuracy of location results by using travel-time based methods. Therefore, many methods that utilize the waveform amplitude have been developed in the past decade. The most representative method is source scanner algorithm (SSA) developed by Kao and Shan (2004, 2007). This method uses the waveform energy in a certain time window to construct an objective function called

ACCEPTED MANUSCRIPT

‘brightness function’ for each grid point in the model. The point which enables the ‘brightness function’ to reach its maximum is assumed to be the event location.

RI PT

Grigoli combined the STA/LTA and ‘brightness function’ together to get a more accurate location result (Grigoli et al., 2014). These methods can be viewed as the combination of ray tracing technique and waveform coherence, therefore, have a high

SC

computation efficiency and noise resistance.

Fink (1999) develops a reverse time imaging technique to back propagate the received

M AN U

waveform at the receivers to focus at its real origin. His method was first demonstrated to be effective in some medical experiments. Compared to reverse time migration in surface seismic exploration, reverse time imaging technique used in source location only has the back propagated wavefield while the migration is

TE D

achieved by cross-correlation of forward propagated wavefield from sources to receivers and the back propagated wavefield from receivers to sources. This method is then applied to locate earthquake by Gajewski (2005), where the point with largest

EP

energy during the back propagation is considered as the earthquake location and the

AC C

corresponding time step is the trigger time. Compared to travel-time based method, this method is more attractive because it follows the dynamic and kinetic characteristics of wave propagation. To perform this method, one has to solve the wave equation to get the back propagated wavefiled. To avoid the large amount of calculations in reverse time technique, Drew (2005) employs the phase shift method instead of finite difference method to get the wavefield at different time and uses a cross-correlation to locate the earthquake. Artman et al. (2010) utilizes the

ACCEPTED MANUSCRIPT

cross-correlation of P wave and S wave energy of back propagated wavefild to estimate the source mechanism. Witten and Artman (2011) develop a noise removal

RI PT

technique for the reverse time imaging to enhance its performance in low S/N measurements. Wang et al., (2013) utilize the Wigner Distribution Function (WDF) to filter the received waveform to suppress the random noise. This technique is similar to

SC

the interferometry imaging technique. To improve the resolution of reverse time imaging, Haldorsen et al. (2012, 2013) uses the deconvolution before imaging.

M AN U

Equation 7 is the acoustic wave equation used in the back propagation. V is the velocity, (x, y, z) is the coordinates of grid point, p is the displacement, D is the received waveform at the receivers, T is the total recording time, t is the time step in reverse time propagation.

1 ∂p( x, y, z, t ) = D( x, y, z = 0, T − t ) 2 v ( x, y , z ) ∂t 2

TE D

∇ 2 p ( x, y , z , t ) +

(7)

Figure 13 shows a 3D numerical example of micro-seismic event location with single

EP

down-hole monitoring data by using reverse time imaging technique. Figure 13 (a)

AC C

shows a four-layered velocity model with one vertical receiver array (red triangular) and a single source (blue dot). The model has 200 grid points in each direction and grid spacing equals to 1 m. The back propagated wavefiled is obtained by 4-order staggered finite difference solver of elastic wave equation and the P and S velocity ratio is a constant of 1.7. Figure 13 (b) shows the reversed imaging in East-North view. The red stars are receivers and the red circle shows the distribution of potential locations. Figure 13 (c) shows the polarization direction of P wave obtained from

ACCEPTED MANUSCRIPT

waveform of three-component receivers. Since the polarization analysis may contain some error due to the noise, here we use Gaussian beam ray to represent the

RI PT

uncertainty. As the distance away from receiver increases, the ray becomes wider, which gives a high uncertainty. Figure 13 (d) shows the East-Depth view of the

reverse time imaging where the virtual source and the real source are displayed

SC

clearly as an energy focus. The virtual source is caused by the 180 degree ambiguity

of the P wave polarization analysis. Figure 13 (e) shows the East-North view of the

M AN U

reverse time imaging combined Gaussian beam ray where the angle uncertainty has been reduced into two directions. Figure 13 (f) shows the North-Depth view. Since the

AC C

EP

TE D

model is symmetrical, the image is same with Figure 13 (d).

Figure 13 Down-hole Micro-seismic event locations using reverse time imaging. (a) a four-layered velocity model with the vertical receiver array (red triangular) and a source (blue dot). Color bar represents the velocity value, which varies from 3000m/s to 4500m/s; (b) the reversed imaging in East-North view; (c) the polarization direction of P wave obtained from waveform of three-component receivers; (d) the East-Depth view of the reverse time imaging; (e) the East-North view of the reverse time imaging combined Gaussian beam ray; (f) the North-Depth view. Colorbars for (b) to (f) show the normalized energy of reversed time images from 0 to 1.

ACCEPTED MANUSCRIPT

The highest energy in the images corresponds to the event location.

Figure 14 shows another 3D numerical example of micro-seismic event location with

RI PT

surface monitoring by using reverse time imaging technique. Figure 14 (a) shows a four-layered velocity model. Figure 14 (b) shows the 11 surface receiver array, and each array contains 11 receivers (blue triangular) and a single source (red star). Figure

SC

14 (c) shows the three component waveform. Figure 14 (d) and (e) show the East North and North-Depth view of reverse time images. It is obvious that the surface

M AN U

monitoring has a better resolution at East-North plane while down-hole monitoring has a better resolution at depth direction.

Here we only show an example of refocusing a point source to its real origin. The influence of different source types especially shear events need to be investigated

TE D

further since the source types will affect the energy distribution in back-propagated wavefield. However, the energy will finally converge at the true location if the velocity model is correct with enough data collected by fully covered receiver

EP

geometry. Artman (2011) use the cross-correlation of P and S wave energy separated

AC C

from the back propagated wavefield as an image condition in reverse time imaging. Different source types, such as the horizontal single point force source, vertical single point force source, 45o single point force and douple couple source are successfully

located and recognized although the energy distribution of final time reversed images are varied.

SC

RI PT

ACCEPTED MANUSCRIPT

M AN U

Figure 14 Surface Micro-seismic event locations using reverse time imaging: (a) a four-layered velocity model; Color bar represents the velocity value, which varies from 3000m/s to 4500m/s. (b) 11 surface receiver array, and each array contains 11 receivers (blue triangular) and a single source (red star); (c) three component waveforms; (d) and (e) are the East - North and North-Depth view of reverse time images, respectively. Colorbars for (d) and (e) show the normalized energy of reversed time images from 0 to 1. The highest energy in the images

TE D

corresponds to the event location.

Compared to travel-time based location algorithms, waveform based methods are able to locate the events without the need of picking events, which greatly enhance our

EP

processing efficiency. In addition, the waveform based methods like reverse time imaging technique employs full waveform data and follows the dynamic and kinetic

AC C

characteristics of wave propagation. In theory, this method is able to refocus the event to its real origin even in strong noise environment. However, the lack of precise velocity model may greatly degrade the accuracy of the location results. In most applications, the velocity model is assumed as horizontal layered model which is commonly acquired from sonic logs or check-shot data. The real velocity of the subsurface may contain strong heterogeneity and anisotropy, which is a great

ACCEPTED MANUSCRIPT

challenge in event location. In addition, with the injection of water, the velocity of the rocks are highly possible to change, therefore, how to consider and model the

interesting issue need to be studied (Zhang et al., 2014).

Relative location method

RI PT

variation of velocity during the injection in real-time event location is another

SC

The location information from methods above is the absolute location and its distribution is diverged, which is very hard to evaluate the geometry. The

M AN U

micro-seismic events are of good similarity if the distance between them is not very far and we can use the cross-correlation method to determine the similarity and the relative arrival times between the events, which can be used for the relative locations between events. The relative location method was proposed by Aster and Scott (1993)

TE D

for earthquake. Since a master event is used to relocate the rest events, the locations are depended on the accuracy of the master event location. Waldhauser and Ellsworth

EP

(2000) addressed a double difference location method which relies on the difference of the arrival times of any two events in the same receiver. The locations obtained

AC C

from this method no longer depend on the accuracy of the master event location.

ACCEPTED MANUSCRIPT

Figure 15. Illustration of the double-difference earthquake relocation algorithm. Solid and open circles represent trial hypocenters that are linked to neighboring events by cross-correlation (solid lines) or catalog (dashed lines) data. For two events, i and j, the initial locations (open circles) and corresponding slowness vectors, s, with respect to two stations, k and l, are shown.

RI PT

Ray paths from the sources to the stations are indicated. Thick arrows (∆x) indicate the relocation vector for events i and j obtained from the full set of equations (9), and dt is the

travel-time difference between the events i and j observed at station k and l, respectively.(Figure

SC

1 in Waldhauser and Ellsworth (2000))

Figure 15 shows the two events i and j and their paths to station k and l. For

M AN U

micro-seismic event i travel to station k, the arrival time
 consists of the original

time   and travel time along the path to station k. It is expressed using ray theory as

a path integral along the ray,


 =   +  



TE D

(8)

Where u is the slowness field and ds is an element of path length. The modeled travel time
_  is usually used a truncated Taylor series expansion (Geiger, 1910) to

EP

linearize. Therefore the difference between the observed and calculated time is,

rki =

∂T ki

∆x li ∑ i ∂x l 3

AC C

=1

+ ∆τ i

(9)

l

For the two nearby events i and j (as shown in Figure 15) are received by station k, the travel time difference between them can be expressed by following equation,

rki − rkj =

∂T i

k ∆x li ∑ i ∂x l 3

=1

l

And rki − rk j = (Tki − Tk j ) obs − (Tki − Tk j ) cal

+ ∆τ i −

∂T j

k ∆x lj ∑ j ∂x l 3

=1

− ∆τ j

(10)

l

consists of the difference of the observed

ACCEPTED MANUSCRIPT

time and the difference of the calculated time for the event pair, which is the origin of the name “double difference”. In traditional location methods, the inaccuracy of

RI PT

arrival time picking and velocity model introduce the main errors for the location. It is clear that the double difference location method can avoid the two factors in some

M AN U

SC

extent.

Figure 16. An example for a micro-seismic events locations by a grid search method (a) and the

TE D

double difference location method (b) (Figure 5 from Zeng et al., 2014)

Here we give an example for micro-seismic events location in a shale gas field of southwestern China by a grid search method (Figure 16a) and the double difference

EP

location method (Figure 16b). The solid curve is the borehole trajectory and points

AC C

with different colors are the events location from different hydraulic fracture stages (as shown in color bar). It is very hard to find any trend from the scattered locations in the Figure 16a. However, two groups of events (labeled by I and II) can be found in Figure 16b, which is improved by the double difference method based on the result in Figure 16a. Group I concentrates around the borehole trajectory corresponding to the hydraulic fracturing treatments and group II is far away from the borehole corresponding to other induced events occurred on the pre-exist fault. Although the

ACCEPTED MANUSCRIPT

double difference location method can be used in the surface monitoring data, we still need to think about the feasibility for the borehole monitoring data due to the azimuth

RI PT

angle and anisotropic in the downhole monitoring.

Moment-tensor inversion of Micro-seismic events

SC

In a general form, the displacement d can be expressed by the convolution of Green’s function (G), moment tensor (M) and source time function (S(t)) as follows, (11)

M AN U

d(x,t) = G*M*S(t),

where x is the location of the sensor and t is the time. The moment tensor can be described by a matrix of nine force couples (Aki and Richards, 2002):

  M =   

  

   

AC C

EP

TE D

And we can also use a figure to describe the elements in the matrix as follows,

Figure 17 Elements of the moment tensor as dipoles in the spatial dimensions and the equivalent forces

The matrix is supposed to be symmetrical and only six elements are independent due

ACCEPTED MANUSCRIPT

to the conservation of angular and linear momentum (Eyre and van der Baan, 2015). The off-diagonal elements are the double couple components (DC) to avoid the

RI PT

rotation. The diagonal elements are force couple for the volumetric changes, which can be used to describe an isotropic source (ISO) when three diagonal elements are

equal (volumetric changed), and also can be a compensated linear vector dipole

SC

(CLVD) when one of the diagonal elements is compensated by others and the

summation of the three diagonal elements is zero (no volumetric changed). If we get

M AN U

the value of the moment tensor matrix, we can get the magnitude of the events, the fracture types and directions which allows us to understand the fracturing behavior

AC C

EP

TE D

and evolving stress field.

Figure 18 Far-field P-wave and SV-wave radiation patterns (red and blue are compressions and dilatations) of double-couple (DC), isotropic (ISO), and CLVD sources, plotted using the

equations of Aki and Richards (2002) and modified from Figure 2 in Eyre and van der Baan

(2015). Arrows at the top schematically show the corresponding force systems generating the source mechanisms (black) and the shear forces (red).

The first-arrival polarity (i.e. the P-wave) or the combination of polarity and S/P

ACCEPTED MANUSCRIPT

amplitude ratio is used in the earthquake for the source mechanism determination, which is the simplest and fastest way. With the good coverage of sensors from

RI PT

different azimuths around the event, we can use the polarity of the P-wave at each sensor to inverse the source radiation pattern, in which the upward first motions are

compressional and the downward first motions are dilatational. The solution is often

SC

double-couple in this method. However, the method requires the monitor geometry

has a very good coverage of sensors and also good SNR of P wave. Therefore, it

M AN U

cannot be applied to the micro-seismic events due to the poor coverage of the sensor system especially for the downhole monitoring and the low SNR of the surface monitoring.

Full waveform inversion method

TE D

It is similar as exploration seismic data processing that the full waveform information is used for the inversion when the first-arrival information is not good for processing.

EP

Because the full waveform includes the scatters from different azimuths, therefore it can be used to inverse the moment tensor even only one azimuth available for the

AC C

monitoring system. According to equation (11), the waveform consists of the contributions from three different parts: the source time function, source mechanism, and propagation path (Green’s function). The Green’s function can be calculated by using different methods such as ray tracing and full-waveform simulations if we have an accurate velocity structure of the model. The 3D velocity model can be obtained from the 3D seismic survey. Unfortunately, the velocity model is usually approximated by the logging data (only P velocity available sometime), which

ACCEPTED MANUSCRIPT

introduces the errors for the full waveform modeling. The source time function is also important issue which is seldom reported by the research in micro-seismic. It is easy

RI PT

to use a combination with simple trigonometric functions to approximate the source time function when the frequency of the event is low (around few Hz), which is common in the earthquake study. However, it would be different for the high

SC

frequency micro-seismic events (more than ten Hz). Different inversion methods can

be used such as weight least-squares method and the quality of the inversion results

M AN U

are tested by the evaluation of the misfit between the synthetic waveforms and the field data.

Li et al. (2011a, 2011b) proposed a full waveform inversion method for source mechanism based on grid searching. By constructing a layered velocity model from

TE D

logs, they used a discrete wavenumber integration method (Bouchon and Coutant, 1994) to calculate all the possible Green’s functions by setting only one element to be 1 and other zeroes each time. And they also assumed the source mechanism to be

EP

double-couple dominated and other mechanisms such as ISO and CLVD were ignored. Then only three independent variables, strike, dip and rake for the solution were

AC C

searched in discrete grids.

Because it is difficult to obtain accurate absolute amplitudes due to site effects in many situations, they normalized the filtered observed and modeled waveforms before comparison. By comparing the synthetic waveforms from all the possible solution with the observed data, they use a cost function to control the search, in which the (1) cross-correlation, (2) the normal 2, (3) the polarities and (4) the S/P

ACCEPTED MANUSCRIPT

amplitude ratios between the normalized modeled and observed waveforms are considered. Energy normalization is used in their method, such that the energy of the

RI PT

normalized wave train within a time window adds to unity. Within different weights for those four parts, the maxima value of the cost function gives a best solution of

source mechanism (see, equation 2 in Li et. al., 2011a).The solution not only includes

SC

the strike, dip and rake, but also includes the possible location of the event. The method was used successfully in both near surface and downhole monitoring. Figure

M AN U

19 shows an example in a gas field in Middle-East for the full waveform source mechanism inversion with spare surface stations (5 stations) (Figure 10 in Li et al., 2011b). The left and right of the figures are the comparisons of P and S waves, respectively, in which the red is modeled data and blue is the observed data. The good

TE D

comparison indicates the method works well even for the poor coverage of the sensors and poor SNR for the surface monitoring data. According to Li et al. (2011b), to reduce the velocity uncertainty, some parameters such as polarities of modeled and

EP

observed compressional waves (“+” and “-” as shown in left column of Figure 19), the arrival time shift of the modeled data, and amplitude ratio between P and S

AC C

waveforms (ratio shown in right top of Figure 19) are needed for matching the observed data. More details about the parameter can be found in Li et al. (2011b).

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

‐

Figure 19. Comparison between the modeled and surface observed waveforms of a micro-seismic event by a source mechanism grid search method. Left and right are for the P and

TE D

S waves. The waveforms are normalized by energy. (Figure 10 in Li et al., 2011b)

Figure 20 shows an example of the source mechanism solutions for different micro-seismic events in the Middle East gas field by Li et al. (2011b).The most of the

EP

events primarily have the normal faulting mechanism, some have the strike-slip

AC C

mechanism, and some have a reverse faulting mechanism.

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 20. Source mechanisms of the 40 events inverted from both the surface and borehole networks. The background color in the map indicates the local change in surface elevation with a maximum difference of about 10 m. Different focal mechanisms are grouped in several colors.

TE D

The events and their focal mechanisms determined by the surface network are plotted in the outer perimeter, while the ones by the borehole network are plotted in the inner ring. (Figure 12a in Li et al., 2011b)

EP

Since the method is based on the double-couple solution assumption, the uncertain

AC C

during the inversion introduced by the anisotropic or in homogeneous medium can be minimized and computation can be simplified. However, the ISO and CLVD solution are not considered, which are also important for understanding the fracturing and evaluating the efficiency of the hydraulic facture operation. Song and Toksoz (2011) determined the full moment tensor solution by a linear inversion method. Although the method is of high effective with very few forward modeling, the accurate location of the event is needed in advance and the P and S waves are hard to distinguish.

ACCEPTED MANUSCRIPT

Imaging in Micro-seismic monitoring technique

To understand the underground structure, the velocity is always a big issue for

RI PT

geoscientist. The velocity of rock will change during the hydraulic fracturing treatments. It is necessary to try to get the velocity information from the

micro-seismic data. Zhang and Thurber (2003) considered the effect of the velocity

SC

model changing on the double difference location method (Waldhauser and Ellsworth,

2001). They use the pseudo-bending ray tracing algorithm (Um and Thurber, 1987) to

M AN U

find the rays and calculate the travel times between events and stations. The model is represented as a regular set of 3D nodes, and the velocity values are interpolated by using the trilinear interpolation method. The hypocentral partial derivatives are calculated from the direction of the ray and the local velocity at the source (Lee and

TE D

Stewart, 1981). In the simultaneous inversion for velocity structure and event locations, velocity anomalies are constrained by seeking a first-order smooth model.

EP

By using their double-difference seismic tomography method, it is possible to get some velocity information around the events. Then we can use the velocity image to

AC C

describe the fractures.

Here, we give an example for velocity imaging from surface micro-seismic monitoring data in a south-west China shale gas field. The P and S velocities in the depth of 2.5 km are imaged in Figure 21 and the events location and the borehole trajectory are shown (white solid line, black are the plug positions). It is clear that almost all the events occur in the low velocity areas which are corresponding to

ACCEPTED MANUSCRIPT

fracturing. It is anticipated that the Vp/Vs plot will help understanding the geological structure. However, the good Vp/Vs plot requires the events with good P and S wave

RI PT

arrival time at the same time. In this field data set, we find that the P wave arrival times of many events are hard to picked which makes the bad quality of the Vp/Vs

M AN U

SC

plot. Therefore, we do not show the Vp/Vs plot in Figure 21.

‐

Figure 21. An example for velocity imaging from surface micro-seismic monitoring data in a

TE D

western China shale gas field: (a) compressional velocity and events location distribution; (b) shear velocity and events location distribution.

From the figure 21, we can find velocity structure is just a smoothed structure which

EP

means the accuracy of the model is not high and the velocity for the area without

AC C

events is inaccuracy. Anisotropic is also a problem for the velocity imaging by the micro-seismic data.

Summary and discussion Recent applications and processing techniques of microseismic monitoring including first arrival picking, absolute and relative location, moment tensor inversion and imaging techniques have been reviewed in this paper. Since microseismic is the only

ACCEPTED MANUSCRIPT

far-field technology that can image the induced fracture geometry within the reservoirs, it is widely applied in hydraulic fracturing monitoring, short-term reservoir

RI PT

characterization, hazard analysis in coal mining and surveillance of gas storage. Accurate and precise processing and interpretation results can help engineers obtain a better understanding of their reservoirs, leading to an optimized production,

SC

stimulation or storage design. However, how to lower the uncertainty and ambiguity

in microseismic monitoring is still a problem although several successful projects

M AN U

demonstrate the potential of this technique.

The current processing technique in industry requires a reliable first-arrival pickings and accurate velocity model. However, an accurate onset picking might be time

TE D

consuming and also very challenging in data with strong random or coherent noise. In recent years, researchers are focusing on techniques that employ not only the arrivals but also the full waveform to locate micro-seismic events by using migration based

EP

methods such as reverse time and diffraction stacking techniques. Unfortunately, a heavy computation cost makes these techniques unrealistic in real-time processing. In

AC C

addition, the migration based methods require a large coverage and dense distribution of sensors. The location results are poorly constrained in a sparse receiver array. Another problem is that the location result suffers a 180 degree ambiguity in single downhole application. The location of injectors is commonly used to judge the direction of the event based on the assumption that the micro-seismic events are tend to occur around the injection intervals. However, it is very difficult to determine the

ACCEPTED MANUSCRIPT

polarization direction of a micro-seismic event in multiple injectors where the depths and azimuth or injection intervals are similar. The most significant problem is the

RI PT

inaccurate velocity model. Sonic logs and VSP provide limited velocity in lateral direction while the surface seismic provides velocity with poor resolution. The lack of

SC

true 3D velocity model significantly affects the reliability of location results.

Although several issues need to be managed, we still consider micro-seismic

M AN U

promising in future. With the development of full waveform inversion (FWI) technique, an accurate velocity model combined with information of anisotropy and attenuation can be obtained. The continuing increase in computer speed and parallel or GPU calculation techniques make it realistic to employ the migration based method.

TE D

A new generation of geophysical measurement called distributed acoustic sensors (DAS) is able to provide complete coverage and dense distribution of sensors both in downhole and surface applications (Sun et al., 2015). 3-C components measuring the

EP

ground velocity plus one additional component measuring the pressure can be applied to reduce the 180 degree ambiguity in single downhole array (Li et al., 2015).

AC C

Therefore, a rapid and automatic processing technique based on massive micro-seismic data might attract the interests of researchers in future. In addition, the microseimisc monitoring technique combined with rock physics, geomechanics and other geophysical measurements such as pressure, temperature as well as injected fluid rates and volumes might greatly constrain the interpretation results with less uncertainty.

ACCEPTED MANUSCRIPT

Acknowledgements This study is supported by an Open Fund of State Key Laboratory of Oil and Gas

RI PT

Reservoir Geology and Exploitation (Southwest Petroleum University, PLN1503), a NSFC (No.41404100), a China Post-doctoral Science Foundation (No.2013M530106)

and The International Postdoctoral Exchange Fellowship Program. The three

SC

anonymous reviewers are deeply appreciated for their critical and constructive

M AN U

comments to this work.

References

Adam Baig and Ted Urbancic (2010). "Microseismic moment tensors: A path to understanding frac growth." The Leading Edge, 29(3), 320-324.

TE D

Aki K., and P. Richards, 2002, Quantitative seismology (2nd ed.): University Science Books. Artman B., I. Podladtchikov and B. Witten, 2010, Source location using time-reverse imaging, Geophysical Prospecting, 2010, 58, 861–873

EP

Allen, R., 1978, Automatic earthquake recognition and timing from single traces, Bulletin of the

AC C

Seismological Society of America, 68, 1521-1532. Aster R. and J. Scott., 1993, Comprehensive characterization of waveform similarity in microearthquake data sets. Bulletin of the Seismological Society of America, 83,1307-1314. Billings S.D., 1994, Simulated annealing for earthquake location, Geophys. J. Int., 118,680-692. Block L.V., C.K. Wood, W.L. Yeck, V.M. King, 2015, Induced seismicity constraints on subsurface geological structure, Paradox Valley, Colorado. Geophysical Journal International 200, 1172-1195.

ACCEPTED MANUSCRIPT

Bouchon M. and O.Coutant, 1994. Calculation of synthetic seismograms in a laterally varying medium by the boundary element-discrete wavenumber method, Bull. seism. Soc. Am., 84(6),

RI PT

1869–1881. Chen Z. and R. R. Stewart, 2006, A multi-window algorithm for real-time automatic detection and picking of P-phases of microseismicevents, CREWES Research Report, Volume 18

SC

Drew J., D. Leslie and P. Armestrong, 2005, Automated Microseismic Event Detection and

Location by Continuous Spatial Mapping. The 2005 SPE Annual Technical Conference and

M AN U

Exhibition, SPE 95513.

Duncan P. M., L. Eisner , 2010, Reservoir characterization using surface microseismic monitoring. Geophysics, 75, 139-146.

Dye B. C., R.H. Jones, J. F. Cowles, O. Barkved and P. G. Folstad, 1999, Microseismic survey of a

TE D

North Sea reservoir. World Oil, 220(3), 74-78.

Eaton D., M. van der Baan, J. Tary, B. Birkelo, N. Spriggs, S. Cuttenand, K. Pike, 2013, Broadband microseismic observations from a Montney hydraulic fracture treatment, northeastern

EP

B.C., Canada,38(3),44-52.

Economides M. J. and K. G. Nolte, 2000, Reservoir Stimulation, John Wiley, West Sussex, U. K.

AC C

Eyre T. and M. van der Baan, 2015, Overview of moment-tensor inversion of micro-seismic events. The leading edge, 34(8),882-888. Finck F., J. H. Kurz, C. U. Grosse, H. Reinhardt, 2003, Advances in Moment Tensor Inversion for Civil Engineering, International Symposium of Non-Destructive Testing in Civil Engineering 2003 Fink M., 1999, Time-reversed acoustics, Scientific American, 281(11), 91–97

ACCEPTED MANUSCRIPT

Gajewski D., E. Tessmer, 2005, Revere modeling for seismic events characterization. Geophysical Journal International, 163,276-284

time only. Bull. St .Louis. Univ, 8, 60~ 71.

RI PT

Geiger L., 1912, Probability method for the determination of earthquake epicenters from arrival

Grigoli F., S. Cesca, O. Amoroso, A. Emolo, A. Zollo, and T. Dahm, 2014, Automated seismic

SC

event location by waveform coherence analysis: Geophysical Journal International, 196(3), 1742−1753.

M AN U

Haldorsen J., M. Milenkovic, N. Brooks, C. Crowell and M.B. Farmani, 2012, A Full-waveform, Migration-based Deconvolution Approach to Locating Micro-seismic Events, 74th EAGE Conference & Exhibition incorporating.

Haldorsen J., N. Brooks and M. Milenkovic, 2013, Migration-based Deconvolution for High

TE D

Resolution Locating of Microseismic Sources from Multi-Component Recordings, Fourth Passive Seismic Workshop Optimizing Development of Unconventional Reservoirs 17-20 March 2013 Hammersley J.M., D.C. Handscomb, 1967, Monte Carlo Methods, Methuen, London.

EP

Healy J. H., W. W. Rubey, D. T. Griggs, C. B. Raleigh, 1968, The Denver earthquakes. Science 161, 1301-1310.

AC C

Hornbach, M.J., DeShon, H.R., Ellsworth, W.L., Stump, B.W., Hayward, C., Frohlich, C., Oldham, H.R., Olson, J.E., Magnani, M.B., Brokaw, C., Luetgert, J.H., 2015, Causal factors for seismicity near Azle, Texas. Nature Commun., 6, 6728. Improta L., L. Valoroso, D. Piccinini, C. Chiarabba, 2015, A detailed analysis of wastewater-induced seismicity in the Val d'Agri oil field (Italy). Geophysical Research Letters, 42, 682-2690.

ACCEPTED MANUSCRIPT

Kalman, R.E., 1960, A new approach to linear filtering and prediction problems: Transaction of the ASME-Journal of Basic Engineering, 33-35.

RI PT

Kamei R., N. Nakata, D. Lumley, 2015, Introduction to microseismic source mechanisms. The Leading Edge, 34(8):876-880

Kanwaldip S. A. and F. U. Dowla, 1997, Wavelet Transform Methods for Phase Identification in

SC

Three-Component Seismograms. Bulletin of the Seismological Society of America, 87(6), 1598-1612

M AN U

Kao H. and S. Shan, 2004, The source-scanning algorithm: mapping the distribution of seismic sources in time and space, Geophys. J. Int., 157,589–594.

Kao H. and S. Shan, 2007, Rapid identification of earthquake rupture plane using source-scanning algorithm, Geophys. J. Int., 168, 1011–1020.

TE D

Langenbruch C. and S. A. Shapiro, 2015, Quantitative analysis of rock stress heterogeneity: Implications for the seismogenesis of fluid-injection-induced seismicity, Geophysics, 80(6), 73-88.

B. Jiang, 2013, A detailed view of the

EP

Lei X., S. Ma, W. Chen, C. Pang, J. Zeng,

injection-induced seismicity in a natural gas reservoir in Zigong, southwestern Sichuan Basin,

AC C

China. Journal of Geophysical Research: Solid Earth, 118, 4296-4311. Leo Eisner, B. J. Hulsey, Peter Duncan, Dana Jurick, 2009, Heigl Werner and William Keller,Comparison of surface and borehole locations of induced seismicity, Geophysical Prospecting, 1, 12. Li J., H. Zhang, H. S. Kuleli and M. N. Toksoz, 2011a, Focal mechanism determination using high-frequency waveform matching and its application to small magnitude induced earthquakes.

ACCEPTED MANUSCRIPT

Geophys. J. Int. ,184 (3), 1261-1274. Li J., H.S. Kuleli, H. Zhang, and M.N. Toksoz, 2011b, Focal mechanism determination of induced

Geophysics, 76(6),87-101.

RI PT

micro-earthquakes in an oil field using full waveforms from shallow and deep seismic networks.

Li Z.H and Mirko, 2015, Enhanced microseismic event localization by reverse time extrapolation.

SC

The 2015 SEG Annual Technical Conference and Exhibition.

Lienert B. R. , Berg E. , Frazer L. N., 1986, Hypocenter: An earthquake location method using

M AN U

centered scaled and adaptively damped least squares. Bull. Seism. Soc. Am, 76(3), 771~783. Maeda, N. (1985). A method for reading and checking phase times in autoprocessing system of seismic wave data, ZisinJishin, 38, 365–379.

Martínez-Garzón P., G. Kwiatek, H. Sone, M. Bohnhoff, G. Dresen, C. Hartline, 2014,

TE D

Spatiotemporal changes, faulting regimes, and source parameters of induced seismicity: A case study from The Geysers geothermal field. Journal of Geophysical Research, 119, 8378-8396. Maxwell S.C., Rutledge, J., Jones, R., and Fehler, M., 2010, Petroleum reservoir characterization

EP

using downhole microseismic monitoring. Geophysics, 75(5), 75A129-75A137. Maxwell S.C., 2011, Microseismic hydraulic fracture imaging: The path toward optimizing shale

AC C

gas production, The Leading Edge, 30, 340-346. Miyazawa, M., Venkataraman, A., Snieder, R., and Payne A., 2008, Analysis of microearthquake data at Cold Lake and its applications to reservoir monitoring. Geophysics 73(3), O15-O21. Nelson G. D., John E Vidale,1990, Earthquake locations by 3-D finite difference travel times. Bull. Seism. Soc. A m , 80(2), 395~ 410. Pearson, C., 1981, The relationship between microseismicity and high pore pressures during

ACCEPTED MANUSCRIPT

hydraulic stimulation experiments in low permeability granitic rocks: J. Geophys. Res., 86, 7855-7864.

RI PT

Prugger A. F., Gendzwill D.J., 1988, Microearthquake location: A nonlinear approach that makes use of a simplex stepping procedure. Bull.Seism.Soc.Am, 78(2), 799~815.

Rutledge, J., Yu, X., Leaney, S., 2015. Microseismic shearing driven by hydraulic-fracture

SC

opening: An interpretation of source-mechanism trends. The Leading Edge 34, 926-934.

Sethian J. and A. Popovici, 1999, 3-D traveltime computation using the fast marching method.

M AN U

Geophysics, 64(2), 516-523.

ílený, J., D. P. Hill, L. Eisner, an d F. H. Cornet, 2009, Non double couple mechanisms of microearthquakes induced by hydraulic fracturing, J. Geophys. Res., 114, B08307, Skoumal, R.J., Brudzinski, M.R., Currie, B.S., 2015a. Earthquakes Induced by Hydraulic

TE D

Fracturing in Poland Township, Ohio. Bulletin of the Seismological Society of America. Skoumal, R.J., Brudzinski, M.R., Currie, B.S., 2015b. Microseismicity induced by deep wastewater injection in Southern Trumbull County, Ohio. Seismological Research Letters.

EP

Sleefe G.E., N.R. Warpinski, and B.P. Engler et al.,1995, The use of broadband microseisms for hydraulic fracture mapping, SPE 26485.

AC C

Sleeman, R., and T. van Eck, 1999, Robust automatic P-phase picking: an on-line implementation in the analysis of broadband seismogram recordings. Phys. Earth Planet. Interiors 113, 265–275. Song F., H.S. Kuleli, M.N. Toksoz, E. Ay, and H. Zhang, 2010, An improved method for hydrofracture-induced microseismic event detection and phase picking. Geophysics, 75(6),47-52.

Song, F., and M. N. Toksöz, 2011, Full-waveform based complete moment tensor inversion and source parameter estimation from downhole microseismic data for hydrofracture

ACCEPTED MANUSCRIPT

monitoring: Geophysics, 76(6), 103–116. Stabile T.A., A. Giocoli, A. Perrone, S. Piscitelli, V, Lapenna, 2014, Fluid injection induced

(southern Italy). Geophysical Research Letters, 41, 5847-5854.

RI PT

seismicity reveals a NE dipping fault in the southeastern sector of the High Agri Valley

Sun J., Zhu T., and Fomel S., 2015, Investigating the possibility of locating microseismic

SC

sources using distributed sensor networks, The 2015 SEG Annual Technical Conference and Exhibition.

M AN U

Um, J., and C. H. Thurber, 1987, A fast algorithm for two-point seismic ray tracing. Bull. Seism. Soc. Am. 77, 972–986.

Waldhauser F., W. L. Ellsworth, 2000, A double - difference earthquake location algorithm: method and applicationto the Northern Hayward Fault, California. Bull.Seism.Soc.Am, 90(6),

TE D

1353~ 1368.

Wang C.L., Cheng J.B., and Yin C. et al., 2013,Microseismic events location of surface and borehole observation with reverse-time focusing using interferometry technique, Chinese J.

EP

Geophys,56(9),584-597.

Warpinski N.R., S.L. Wolhart, C.A. Wright, 2001, Analysis and prediction of microseismicity

AC C

induced by hydraulic fracturing, SPE 71649. Wei, S., Avouac, J.-P., Hudnut, K.W., Donnellan, A., Parker, J.W., Graves, R.W., Helmberger, D., Fielding, E., Liu, Z., Cappa, F., 2015. The 2012 Brawley swarm triggered by injection-induced aseismic slip. Earth and Planetary Science Letters 422, 115-125. Will, R., Smith, V., Leetaru, H.E., Freiburg, J.T., Lee, D.W., 2014. Microseismic monitoring, event occurrence, and the relationship to subsurface geology. Energy Procedia 63, 4424-4436.

ACCEPTED MANUSCRIPT

Witten B, and B Artman, 2011, Signal-to-noise estimates of time-reverse images, Geophysics, 76(2): MA1–MA10

RI PT

Zeng, X., H. Zhang, X. Zhang, H. Wang, Y. Zhang, Q. Liu, 2014, Surface microseismic monitoring of hydraulic fracturing of a shale gas reservoir using short period and broadband seismic sensors, Seismological Research Letters, 85(3), 668-677.

SC

Zhang H., C. Thruber, C. Rowe, 2003, Automatic P-Wave Arrival Detection and Picking with

Multiscale Wavelet Analysis for Single-Component Recording. Bulletin of the Seismological

M AN U

Society of America, 95(5),1904-1912.

Zhang H., C. Thurber, 2003, Double-difference tomography: The method and its application to the Hayward fault, California. Bulletin of the Seismological Society of America, 93 (5), 1875-1889 Zhang X.L., Zhang F., Li X. Y., et al., 2014, The influence of fracturing process on microseismic

TE D

propagation, Geophysical Prospecting, 62,797-805.

Zhao A H, Ding Z F, Sun W G, 2008, Calculation of focal loci for eqrthquake location in complex media, Chinese J.Geophys, 51(4),845-854

EP

Zhou J C , Zhao A H., 2012, An intersection method for locating earthquakes in 3-D complex velocity models, Chinese J. Geophys, 55(10). 3347-3354

AC C

Zollo A, Capuano P, Virieux J (2001), Precise, absolute earthquake location under Somma-Vesuvius volcano using a new 3D velocity model, Geophys. J. Int., 146,313-331.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT




Briefly review of characteristics of micro-seismic events, event onset time auto-picking, event location and relocation as well as source mechanism analysis. The uncertainties and issues of current processing workflow are addressed




Finally potential improvement and application are discussed.

AC C

EP

TE D

M AN U

SC

RI PT