Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir

- TREBALL FI DE CARRERA Títol Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir Autor/a Marc Gerona Sendino Tutor/a Alan F. Baird Lluis Pujades Beneit Departament Geophysics group Departament d’Enginyeria del Terreny, Cartogràfica i Geofísica Intensificació Data Maig de 2014 Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir I My most sincere thanks To my supervisor, Alan F. Baird, for giving me the opportunity to develop this research project in the geophysics group at the University of Bristol. For his support, patience and for all the knowledge he has taught me during this year. To Lluis Pujades, for helping me whenever I need him and for all the interest he has always shown me. Finally, to my family and friends. And especially to Laura, for giving me all her support from the distance. Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir III Abstract The exploitation of unconventional reservoirs has become an important component for the world’s gas and oil supplies. Hydraulic fracture stimulations are used to generate fractures in the rock, improving the permeability of the reservoirs and increasing the productivity of the wells. Induced fractures generate microseismic events which can be used to monitor the evolution of the fracture network. The first step of a microseismic monitoring is to detect the events using an effective automated method. Here, we use a continuous passive microseismic dataset, recorded during hydraulic fracture stimulations in a shale oil reservoir, to develop a new automated microseismic event detection method. The microseismic acquisition company also provided event sorted data they picked which is used to measure the effectiveness of the new method. Some data gaps identified in the event data provided by the acquisition company have been filled by the new method. Additionally, the new method has detected a large number of events during the first day of monitoring, which were not included in the event data provided. Further analysis has been done using the event data provided by the acquisition company. The b-values and DC-values are estimated to provide an insight into the mode of failure during the hydraulic fracture stimulation and the complexity of the induced fractured network. A simpler fracture network is inferred for stages 2, 3 and 4 with events likely clustered along a planar feature. In contrast, a more complex fracture network is deduced for stage 1 with events may randomly distributed throughout a 3D volume. Many reasons have been contemplated to explain this difference but unfortunately, more information is needed in order to be more precise and confident. Finally, shear wave splitting of event data provided is analysed to estimate seismic anisotropy which could provide insight into the natural and induced fractures. It has been inferred that the anisotropy of the medium is likely dominated by near-vertical fractures oriented in the direction of maximum stress (NE-SW). Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir V Table of contents Thanks ........................................................................................................................................... I Abstract ...................................................................................................................................... III Table of contents .......................................................................................................................... V List of figures ............................................................................................................................ VII List of tables ................................................................................................................................. X 1. Introduction ................................................................................................................ 1 1.1. Objectives ......................................................................................................................... 2 1.2. Outline .............................................................................................................................. 3 2. Microseismic monitoring context .............................................................................. 4 2.1. Geological context ........................................................................................................... 4 2.2. Well completion and orientation ...................................................................................... 5 2.3. Locations and settings of the wells and the geophones .................................................... 8 3. Microseismic data ..................................................................................................... 12 3.1. Continuous data .............................................................................................................. 12 3.2. Event data ....................................................................................................................... 14 4. Automated microseismic event detection method ................................................. 17 4.1. STA/LTA algorithm ....................................................................................................... 17 4.1.1. STA/LTA parameters ............................................................................................ 19 4.2. STA/LTA trigger algorithm ........................................................................................... 23 4.2.1. STA/LTA trigger parameters ................................................................................ 23 4.3. Cross-correlation ............................................................................................................ 28 4.4. Results of the new event detection method .................................................................... 35 5. Magnitude-frequency and spatial distributions .................................................... 38 Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir VI 5.1. Magnitude-frequency distribution .................................................................................. 38 5.1.1. Estimation of b-value ........................................................................................... 39 5.2. Spatial distribution ......................................................................................................... 43 5.2.1. Calculation of DC ................................................................................................. 43 5.3. Interpretation of the magnitude-frequency and spatial distributions ............................... 49 6. Shear wave splitting ................................................................................................. 51 6.1. Shear wave splitting measurements ............................................................................... 54 6.2. Interpretation of the results ............................................................................................ 56 7. Conclusions and future research ............................................................................. 59 7.1. Conclusions .................................................................................................................... 59 7.2. Future research ............................................................................................................... 60 8. References .................................................................................................................. 61 Appendix ..........................................................................................................................i Appendix A: Code in Matlab of the automated microseismic event detection method ......... ii Appendix B: Code in Matlab created by James Verdon in 2012 to estimate b-values ....... xxi Appendix C: Code in Matlab created by James Verdon in 2012 to estimate DC-values .. xxvi Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir VII List of figures Figure 2.1. Extension of the Bakken Formation and the Williston Basin in Canada. ................. 4 Figure 2.2. Schematic general chart showing Bakken Formation members. .............................. 5 Figure 2.3. Schematic 3D section of the fracture stimulation in the Bakken formation. ............ 6 Figure 2.4. Schema that shows the elements of the system used to stimulate the microseismic monitored well. ............................................................................................................................ 6 Figure 2.5. World stress map zoom of the zone of the microseismic monitored well. ............... 7 Figure 2.6. Schematic diagram shows the different fracture networks produced by stimulation in two different oriented wells. .................................................................................................... 7 Figure 2.7. P and S waves velocity model sent by the acquisition company. ............................. 8 Figure 2.8. Location of the stimulation zone. .............................................................................. 9 Figure 2.9. Perspective view showing location of the well heads and dynamite shots. .............. 9 Figure 2.10. Schematic diagrams showing the distribution of the elements of the microseismic monitoring from the top and side views. .................................................................................... 10 Figure 3.1. Plot of the three traces of every geophone, green and red for horizontal traces and blue for vertical trace, for a specific data file. ............................................................................ 12 Figure 3.2. A zoom of the graph in the Figure 3.1 between the seconds 4 and 4.5. .................. 13 Figure 3.3. Plot of the sum of the three traces of every geophone for the data File in the Figure 3.1. .......................................................................................................... 13 Figure 3.4. Sum of the three traces of every geophone during an event picked by the company. .............................................................................................................. 14 Figure 3.5. 3D Plot with the string of geophones, the stimulation well and the locations of the events picked by the acquisition company. ................................................................................ 15 Figure 3.6. Scattered graph showing the magnitude versus distance from the receiver for all the events picked by the acquisition company. ................................................................................ 15 Figure 3.7. Histogram of the events picked by the company with 1 minute bins and different colour for every stage. ................................................................................................................ 16 Figure 4.1. The different steps of the STA/LTA algorithm around a recorded event. .............. 18 Figure 4.2. Plot of the sum of the three traces of every geophone of one of the files that have been used to test the behaviour of the STA/LTA algorithm. ..................................................... 21 Figure 4.3. STA/LTA ratio functions obtained applying the STA/LTA algorithm, using the different approaches and sets of parameters in the Tables 4.1 and 4.2, to the trace of the Geophone 6 of the file plotted in the Figure 4.2. ....................................................................... 22 Figure 4.4. Plot of the sum of the three traces of every geophone of the first file we have used to show how works the triggering algorithm. ............................................................................. 23 Figure 4.5. Plot of the geophone 8 trace from the file plotted in the Figure 4.4. ...................... 24 Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 4 2. MICROSEISMIC MONITORING CONTEXT 2.1 Geological context The microseismic dataset was acquired during a multi-stage hydraulic fracture stimulation in Southeast Saskatchewan. The oil was produced from the Bakken Formation, which is restricted to the sub-surface in the Williston Basin in southern Saskatchewan and south-western Manitoba in Canada; and in western North Dakota and north-eastern Montana in the United States (Figure 2.1). Figure 2.1. Extension of the Bakken Formation and the Williston Basin in Canada [SK (Saskatchewan), MB (Manitoba)] and the United States [MT (Montana), ND (North Dakota) and SD (South Dakota)] and location of the multi-stage hydraulic fracture stimulation (http://www.undeerc.org/Bakken/bakkenformation.aspx). The Bakken Formation is a part of a vast interval of Late Devonian and early Mississippian black shale formations. It is composed of two hemipelagic mudstone members (upper and Multi-stage hydraulic fracture stimulation Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 5 lower) separated by a shallow marine, grey mudstone/sandstone middle member (Smith and Bustin, 1998) (Figure 2.2). Porosity and permeability within the middle member are generally very low. Porosity averages 5% while permeability averages 0.04x10-3 millidarcies. As the organic content of the middle member is generally very low, it serves as trap for oil derived from the upper and lower members, considered source rocks of the formation (Pitman, 2001). The middle member is considered a tight reservoir due to its low permeability and porosity so it needs to be stimulated in order to become productive. . Figure 2.2. Schematic general chart showing Bakken Formation members (http://www.undeerc.org/Bakken/bakkenformation.aspx). 2.2 Well completion and orientation Generally, stimulation wells consist of a vertical section in which the desired depth is reached and a horizontal section where the hydraulic fracturing is done. In the Bakken Formation the horizontal section of the wells are completed in the upper middle Bakken member in order to increase the effective conductivity of oil from the upper and lower Bakken shale (Figure 2.3). There is significant variability in the well completions across the Williston Basin. In this case, the microseismic monitoring was developed in a lateral horizontal well completed with a system called mechanical isolators. The system comprises ported sleeves installed between isolation packers on a single liner string (Figure 2.4). The ported sleeves activate the packers progressively during the hydraulic fracturing. Packers isolate the horizontal wellbore into Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 6 different stages to stimulate different zones of the formation. The well is cleaned out by flow back to the surface which returns fluid and solid particles. Figure 2.3. Schematic 3D section of the fracture stimulation in the Bakken formation. Figure 2.4. Schema that shows the elements of the system used to stimulate the microseismic monitored well (http://www.undeerc.org/Bakken/bakkenformation.aspx). The region of the Bakken formation is dominated by northeast-southwest maximum stress orientations as can be seen in the Figure 2.5. In order to take advantage of induced fracture propagation in the direction of maximum stress, most wells are positioned in a north-south or northwest-southeast orientation (Figure 2.6). The orientation of the stimulation well monitored is approximately north-south. Bakken Formation Induced fractures Stimulation well Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 7 Figure 2.5. World stress map zoom of the zone of the microseismic monitored well (Heidbach et al. The world Stress Map database release 2008). Figure 2.6. Schematic diagram shows the different fracture networks produced by stimulation in two different oriented wells (Wright et al. 2013). Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 8 1350 1400 1450 1500 1550 1600 1650 1700 1750 1800 1600 2600 3600 4600 5600 Depth (m) Velocity (m/s) P-velocity (m/s) S-velocity (m/s) 2.3 Locations and settings of the wells and the geophones The company sent valuable information related with the microseismic monitoring, including a microseismic quality control report, a velocity model and excel sheets with information of the wells and the geophones. Analysing the velocity model of the acquisition company (Figure 2.7) it can be seen that the velocity of the S and P is stable unless between 1450 and 1550 meters where the P and S velocities have a significant decrease. Figure 2.7. P and S waves velocity model sent by the acquisition company. The stimulation well monitored is located in southeast Saskatchewan (Canada), 120 kilometers to the south-east from Regina, very close to an industrial town called Creelman (Figure 2.8). The microseismic monitoring is developed in an almost vertical well located near the zone of the stimulation. A string of twelve 3-component geophones was installed downhole in the observer well, at the Bakken Formation depth. The orientation of the geophones was determined by three shallow dynamite charges (1/4 Kg at 12m) similar to those used in conventional land seismic acquisition. Locations of the stimulation well, the observer well and the dynamite shots are shown in the Table 2.1 and the Figure 2.9. Two schematic diagrams showing the elements distribution of the microseismic monitoring are represented in the Figure 2.10. Finally, all the settings of the twelve 3-component geophones string are represented in the Table 2.2. Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 9 Observer wellhead Stimulation wellhead Shallow charges Figure 2.8. Location of the stimulation zone. Wells Location relative to Observer wellhead Elevation (m) S(-)/N(+) (m) W(-)/E(+) (m) Stimulation wellhead 1311 -87 613 Dynamite trench 1 91 183 617 Dynamite trench 2 -3 182 617 Dynamite trench 3 -99 183 617 Observer wellhead Location (Geographical coordinates) 617 Latitude Longitude 49° 43' 48.4'' N 103° 18' 51.2'' W Table 2.1. Location of the stimulation wellhead and the dynamite explosions relative to the geographical coordinates of the Observer wellhead. Figure 2.9. Perspective view showing location of the well heads and dynamite shots. Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 10 SIDE VIEW 12, 3-Component Geophones Bakken reservoir Stimulation well (final section) Observer well section ~ 10.4 m TOP VIEW 90 m Stimulation well (final section) Observer wellhead Shallow charges . Figure 2.10. Schematic diagrams showing the distribution of the elements of the microseismic monitoring from the top and side views. Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 11 Table 2.2. Settings of the string of twelve 3-component geophones that was installed downhole in the monitor well, at the Bakken Formation profundity. *SW (Stimulation wellhead) **OW (Observer wellhead) . Orientation of the geophones (°) 66.57 70.79 89.11 74.91 53.49 57.42 78.54 63.24 -205 -91.3 76.27 -76.7 Direction of inclination (°) 39.3 39.2 39 38.9 38.7 38.6 38.4 38.3 38.1 38 37.9 37.8 Inclination (°) 0.7 1.1 1.5 1.8 2.2 2.6 3 3.4 3.8 4.2 3.9 3.3 Location respect the OW** S(-)/N(+) (m) 26.65 26.91 27.17 27.44 27.7 27.97 28.23 28.5 28.76 29.02 29.12 29.17 W(-)/E(+) (m) -21.78 -21.88 -21.98 -22.07 -22.17 -22.27 -22.36 -22.46 -22.55 -22.65 -22.65 -22.63 Depth respect the SW* (m) 1450.1 1460.5 1470.8 1481.2 1491.6 1502.1 1512.5 1522.9 1533.3 1543.7 1554.1 1564.5 Space (m) 10 10 10 10 10 11 10 10 10 10 10 - Geophones 1 2 3 4 5 6 7 8 9 10 11 12 Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 12 0 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 11 12 Time (seconds) Geophones File000067.dat [ 02-oct-2007 / 12:01:47 - 12:01:57 (Local time) ] 3. Microseismic information The microseismic monitoring was active during three non consecutive days: 2nd, 3rd and 5th of October in 2007. As a result, the company obtained a large continuous microseismic dataset which was stored in SEG2 format. The acquisition company went through the continuous data detecting events automatically for different stages, creating one file for every event detected. A Matlab code called “seg2_read” created by James Wookey in 2008 is used to read the files that contain the data. The main function of this code is reading the numeric data of every trace. In addition, this code shows you metadata of the file as the acquisition date, the acquisition time, the company, the client, the instruments and the units used, the number of traces in the file and particular information of each trace, for instance the sample interval. 3.1 Continuous data Using the Matlab code “seg2_read” to read some files of the continuous data important information is obtained. As expected, every file has 36 traces due to the company used twelve 3component geophones, so every geophone has three traces, two horizontal and one vertical. Every trace has the same sample interval, 0.00025 seconds, and the same number of points, 40.000. Therefore every file contains 10 seconds of microseismic data for every geophone. The different geophone traces of a continuous file have been plotted overlaid in the Figure 3.1. Figure 3.1. Plot of the three traces of every geophone, green and red for horizontal traces and blue for vertical trace, for a specific data file. Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 13 4 4.05 4.1 4.15 4.2 4.25 4.3 4.35 4.4 4.45 4.5 1 2 3 4 5 6 7 8 9 10 11 12 Time (seconds) Geophones File000067.dat [ 02-oct-2007 / 12:01:47 - 12:01:57 (Local time) ] 0 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 11 12 Time (seconds) Geophones File000067.dat [ 02-oct-2007 / 12:01:47 - 12:01:57 (Local time) ] Despite the data is quite clean and so the signal-noise ratio is high, a bandpass filter of [75Hz 350 Hz] is applied. It can be seen that the graph is quite confusing so a zoom is needed (Figure 3.2) to see the three traces of every geophone. In order to avoid this problem, the sum of the three traces of every geophone is going to be plotted henceforth to get a better visualization of the data (Figure 3.3). Figure 3.2. A zoom of the graph in the Figure 3.1 between the seconds 4 and 4.5. Figure 3.3. Plot of the sum of the three traces of every geophone for the data file in the Figure 3.1. Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 20 Approaches to apply the algorithm to the data A Sum the three components of the geophone and get the STA/LTA ratio waveform B Get the STA/LTA ratio waveforms of every component of the geophone and sum them C Get the STA/LTA ratio waveforms of every component of the geophone and multiply them Table 4.1. Three different approaches to apply the STA/LTA algorithm to the data. The number of data samples of the short-term average windows should be 2-3 times the dominant period of the signal, and the number of data samples of the long-term average windows should be 10 times or more the number of data samples of short-term average windows to produce good STA/LTA results (Akram and Eaton, 2012). Five sets of parameters are proposed in Table 4.2, where only sets 2, 3 and 4 have been chosen based on the Akram and Eaton statement. The behaviour of the algorithm is tested applying it to many files of the data in the three different approaches with the five different sets of parameters. Some graphs are plotted to show the results obtained for a specific geophone of a file (Figure 4.2 and 4.3). Sets of parameters 1 2 3 4 5 Nº samples Short Time Average windows (S) 100 200 300 400 500 Nº samples Long Time Average windows (L) 500 2500 3000 4500 10000 Table 4.2. Five sets of STA/LTA algorithm parameters proposed for test the behaviour of the algorithm. Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 21 0 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 11 12 File000021.dat [ 02-oct-2007 / 11:54:07 - 11:54:17 (Local time) ] Time (seconds) Geophones Figure 4.2. Plot of the sum of the three traces of every geophone of one of the files that have been used to test the behaviour of the STA/LTA algorithm. Analysing the results (Figure 4.3) it can be seen that the sets of parameters 1 and 5 produce poor results. The STA/LTA waveform obtained with the first set is quite noisy in all the approaches to apply the algorithm which makes difficult to pick the event correctly. The STA/LTA waveforms obtained with the fifth set are quite good to pick the event. However, the algorithm is not taking into account one quarter of the trace due to the high size of the long time average window, so a lot of events located in this part of the trace could be missed. As expected, the sets 2, 3 and 4 which follow the guidelines proposed by Akram and Eaton produce good STA/LTA results. Calculating the STA/LTA ratio functions of every component and multiplying them provides a better contrast between the event and the rest of the trace. The set 2 produces the highest STA/LTA ratio values in the zone of the event and is the set with the lowest size of the long time average window, which means that a little part of the file is going to be missed. As picking the events is easier when the STA/LTA ratio values are much higher in the zone of the events respect the rest of the trace, the algorithm is going to be applied using the second set of parameters to produce the STA/LTA ratio functions of every component and multiplying them together. Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 22 0 1 2 3 4 5 6 7 8 9 10 0 2 4 STA/LTA ratios from Geophone 6 trace ( File000021.dat ) / A = Way the algorithm is applied and (number) = Set of parameters used A (1) 0 1 2 3 4 5 6 7 8 9 10 0 5 10 A (2) 0 1 2 3 4 5 6 7 8 9 10 0 5 10 A (3) 0 1 2 3 4 5 6 7 8 9 10 0 5 10 A (4) 0 1 2 3 4 5 6 7 8 9 10 0 5 10 A (5) Time (seconds) STA / LTA ratios 0 1 2 3 4 5 6 7 8 9 10 0 5 10 STA/LTA ratios from Geophone 6 trace ( File000021.dat ) / B = Way the algorithm is applied and (number) = Set of parameters used B (1) 0 1 2 3 4 5 6 7 8 9 10 0 10 20 B (2) 0 1 2 3 4 5 6 7 8 9 10 0 10 20 B (3) 0 1 2 3 4 5 6 7 8 9 10 0 10 20 B (4) 0 1 2 3 4 5 6 7 8 9 10 0 10 20 30 B (5) Time (seconds) STA / LTA ratios 0 1 2 3 4 5 6 7 8 9 10 0 10 20 30 STA/LTA ratios from Geophone 6 trace ( File000021.dat ) / C = Way the algorithm is applied and (number) = Set of parameters used C (1) 0 1 2 3 4 5 6 7 8 9 10 0 100 200 300 C (2) 0 1 2 3 4 5 6 7 8 9 10 0 100 200 C (3) 0 1 2 3 4 5 6 7 8 9 10 0 100 200 C (4) 0 1 2 3 4 5 6 7 8 9 10 0 500 1000 C (5) Time (seconds) STA / LTA ratios Figure 4.3. STA/LTA ratio functions obtained applying the STA/LTA algorithm, using the different approaches and sets of parameters in the Tables 4.1 and 4.2, to the trace of the Geophone 6 of the file plotted in the Figure 4.2. Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 23 0 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 11 12 File000014.dat [ 02-oct-2007 / 11:52:57 - 11:53:07 (Local time) ] Time (seconds) Geophones 4.2 STA/LTA trigger algorithm A STA/LTA trigger algorithm is used to detect the events of the data. The STA/LTA trigger algorithm has four parameters: high threshold, low threshold, minimum triggering gap and maximum triggering gap. The high threshold determines the limit, in which the trigger is on, that means a possible start of an event is detected, and the low threshold determines the limit in which the trigger is off, that means a possible end of an event is detected. The minimum triggering gap is the minimum interval of time allowed between two different triggered events and the maximum triggering gap is the maximum interval of time between the start and the end of the same event. 4.2.1 STA/LTA trigger parameters The values of the thresholds and the gaps have to be determined by trial and error. A lot of files of the data were studied in order to determine the thresholds, though only the results of two of them are shown. The first file has an event between the seconds 8 and 10 (Figure 4.4). The trigger algorithm is applied using the STA/LTA ratio function of the geophone 8, whose signal trace can be seen in the Figure 4.5, where a small signal at the beginning and the event at the end are appreciated. The STA/LTA ratio function has a small increase at the beginning and a large increase at the end, representing the small signal and the event respectively (Figure 4.6). Giving a value of 200 to the high threshold and a value of 0.5 to the low threshold, the STA/LTA trigger algorithm detects the event, as can be seen in the Figures 4.6 and 4.7. It is important to say that this process is made for every geophone of the file, but it is easier and clearly to show only the results from one of them. Figure 4.4. Plot of the sum of the three traces of every geophone of the first file we have used to show how works the triggering algorithm. Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 24 0 1 2 3 4 5 6 7 8 9 10 0 200 400 600 800 1000 1200 1400 1600 Time (seconds) STA / LTA ratio STA / LTA ratio from Geophone 8 trace (File000014.dat) with high and low thresholds STA / LTA waveform High threshold Low threshold 0 1 2 3 4 5 6 7 8 9 10 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 Geophone 8 trace (File000014.dat) [ 02-oct-2007 / 11:52:57 - 11:53:07 (Local time) ] Time (seconds) Amplitude Figure 4.5. Plot of the geophone 8 trace from the file plotted in the Figure 4.4. Figure 4.6. STA/LTA ratio function for the geophone 8 trace plotted in the Figure 4.5 and chosen thresholds. Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 25 0 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 11 12 Time (seconds) Geophones File000025.dat [ 02-oct-2007 / 11:54:47 - 11:54:57 (Local time) ] 0 1 2 3 4 5 6 7 8 9 10 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 Event picked in the Geophone 8 trace (File000014.dat) [ 02-oct-2007 / 11:52:57 - 11:53:07 (Local time) ] Time (seconds) Amplitude Figure 4.7. Plot of the geophone 8 trace of the Figure 4.5 with a mark where the trigger algorithm is on (red) and off (yellow). Analysing the results, looks like the triggering is working quite well because it has picked the event we had observed. However, if the same process is repeated for the geophone 7 of the second file, which has an event between the seconds 5 and 7 (Figures 4.8 and 4.9), there is a problem because the amplitude of the signal in this case is much lower than the signal of the first file, as can be seen comparing the Figures 4.5 and 4.9. Due to the low amplitude of the signal, the values of the STA/LTA ratio function are much lower so the high threshold of 200 is too high (Figure 4.10). Figure 4.8. Plot of the sum of the three traces of every geophone of the second file we have used to show how works the triggering algorithm. Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 26 0 1 2 3 4 5 6 7 8 9 10 0 50 100 150 200 250 STA / LTA ratio from Geophone 7 trace (File000025.dat) with high and low thresholds Time (seconds) STA / LTA ratio STA / LTA wavefrom High threshold Low threshold 0 1 2 3 4 5 6 7 8 9 10 -8 -6 -4 -2 0 2 4 6 8x 10-3 Time (seconds) Amplitude Geophone 7 trace (File000025.dat) [ 02-oct-2007 / 11:54:47 - 11:54:57 (Local time) ] Figure 4.9. Plot of the sum of the geophone 7 traces from the file plotted in the Figure 4.8. Figure 4.10. STA/LTA ratio function for the geophone 7 trace plotted in the Figure 4.9 and chosen thresholds. Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 27 0 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 11 12 Events picked in the File000014.dat [ 02-oct-2007 / 11:52:57 - 11:53:07 (Local time) ] with high threshold = 50 Time (seconds) Geophones 0 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 11 12 Events picked in the File000025.dat [ 02-oct-2007 - 11:54:47 - 11:54:57 (Local time) ] with high threshold = 50 Time (seconds) Geophones Giving a value of 50 to the high threshold the STA/LTA trigger algorithm should pick the event. The results of applying the trigger with a high threshold of 50 to both files are plotted in the Figures 4.11 and 4.12. Although the trigger is detecting both events, it is giving false detections as well. That means that the trigger is also picking false events or small signals that do not have continuity through the other geophones. Even for the first high threshold of 200 the trigger is giving a false detection in the first file studied, as can be seen in the Figure 4.13. Figure 4.11. Events picked with a high threshold of 50 in the first file studied ‘File000014.dat’. Figure 4.12. Events picked with a high threshold of 50 in the second file studied ‘File000025.dat’. Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 28 0 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 11 12 Events picked in the File000014.dat [ 02-oct-2007 11:52:57 - 11:53:07 (Local time) ] with high threshold = 200 Time (seconds) Geophones Figure 4.13. Events picked with a high threshold of 200 in the first file studied ‘File000014.dat’. When an automated code to pick events from a data is being developed, it is difficult to find a balance in the triggering parameters in order to pick only real events. Therefore, you try to pick as much real events as you can, minimising the false detections. As we have seen, with higher values of the high threshold the trigger has less false detections but miss events with low amplitudes. We want to pick the maximum number of real events, so we must figure out how to reduce the number of false events picked keeping the value of the high threshold in 50. The values for the gaps have been determined taking into account the duration of the events and the time between them in all the files studied. The maximum triggering gap has a value of 5000 samples (1.25 seconds) and the minimum triggering gap has a value of 7500 samples (1.875 seconds). 4.3 Cross-correlation It has been found that even for traces with a low signal-to-noise ratio, the waveform of the STA/LTA ratio is similar for all the traces irrespective of their individual signal-to-noise ratios (Forghani-Arani et al. 2013). Therefore, the similarity of STA/LTA functions may be used to track the picked event across geophones in an array. The STA/LTA ratio functions of all the geophones of the first file analysed in the triggering are plotted in the Figure 4.14. As expected, in the zone of the event the STA/LTA ratio functions are similar for almost all the traces, whereas in the zone of the small signals picked as events by Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 29 0 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10 11 12 Time (seconds) Geophones STA/LTA ratio waveforms of the traces from the 'FILE000014.dat' file the triggering threshold, the peak of the function is quite small and looks like there are only similar peaks in two or three additional geophones. Cross-correlation, which is a measure of similarity of two different waveforms as a function of a time-lag applied to one of them, is used in order to quantify the similarity or dissimilarity amongst STA/LTA ratio functions at a particular time (Forghani-Arani et al. 2013). A Matlab code which use cross correlation to declare if the events picked by the trigger are false or real has been developed. The STA/LTA ratio functions of the geophones 3 and 8 in the Figure 4.14 are going to be used to show an example of every step of the Matlab code. The events picked by the triggering threshold in the geophone 8 are going to be declared real or otherwise false. We refer the event picked between the seconds 2 and 3 as event 1 and the event picked between the seconds 8 and 9 as event 2. Two different windows are defined around the event. A small window, calculated resting 1000 samples (0.25 seconds) to the start of the event and adding 1000 samples (0.25 seconds) to the end of the event. A big window, calculated in the same way but resting 2500 samples (0.625 seconds) and adding 2500 samples (0.625 seconds) to the start and the end of the event respectively. The small window is used to cut the STA/LTA ratio function of the geophone where the event has been picked and the big window to cut the functions of the remaining geophones (Figure 4.15). Figure 4.14. STA/LTA ratio function for every geophone trace in the first file studied ‘File000014.dat’. Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 36 09:30 10:30 11:30 12:30 13:30 08:30 09:30 0 1 2 3 4 5 6 7 8 Local time (03-oct-2007) Local time (05-oct-2007) Number of events (1 minute bins) Histogram of the events detected with the new automated microseismic event detection method in the range of time of the events detected by the acquisition company Events - Stage 1 Events - Stage 2 Events - Stage 3 Events - Stage 4 It is important to compare these results with those obtained by the acquisition company to ensure the good performance of the new method. Using the same time interval of the histogram of the events detected by the acquisition company, a new histogram showing the microseismic events detected by the new method can be seen in the Figure 4.25. Comparing both histograms (Figures 3.6 and 4.25) can be seen that their general shape is quite similar. Furthermore, the new method has detected more events in all the stages as can be seen in the Table 4.3 and in the Figure 4.26. The new method has detected the double number of events in the stages 1 and 3, five times more events in the stage 4 and roughly the same events in the stage 2. The difference in the stages 1 and 3 is largely due to data gaps in the events provided by the acquisition company between 9:30-10:00 and 11:30-12:00 which have been filled by the new method. Figure 4.25. Histogram of the events picked by the new automated event detecting code in the company events interval time with 1 minute bins and different colour for every stage. Acquisition company New detection method Events detected - Stage 1 55 93 Events detected - Stage 2 37 44 Events detected - Stage 3 98 187 Events detected - Stage 4 18 89 Table 4.3. Differences between number of events detected by the new method and by the acquisition company. Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 37 09:30 10:15 11:05 11:55 12:45 09:0508:40 0 50 100 150 200 250 300 350 400 450 Local time (03-oct-2007) Local time (05-oct-2007) Cumulative number of events detected Acquisition company events New automated method events Figure 4.26. Cumulative number of events picked by the acquisition company and by the new automated method during the same interval of time. Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 38 0 1 2 3 4 5 6 7 10-1 100 101 102 103 104 105 106 Moment magnitude Nº events (logarithmic scale) b-value estimation for an event population (East Canada) 5. Magnitude-frequency and spatial distributions A population of microseismic events can be described by its magnitude-frequency and spatial distributions, which can provide insight into the mode of failure during the hydraulic fracture stimulation and the complexity of the induced fractured network. As the location and the moment magnitude of the events is needed for this analysis, the next step would be to estimate them for the events picked by the new method. However, we do not have enough time to do the estimation and the analysis, therefore the magnitude and locations provided by the acquisition company are used. 5.1 Magnitude-frequency distribution The magnitude-frequency distribution for an event population is described by the well-known Gutenberg-Richter relationship: (5.1) Where NM is the number of events with magnitude greater than M, and a and b are constants to be determined. The b-value is the gradient of the magnitude distribution. Natural earthquake populations generally have a b-value around 1.0, as can be seen in the example of the Figure 5.1, where the b value for an event population located in eastern Canada has been estimated. Figure 5.1. Approximation of the b value of an event population located in eastern Canada. 10 log M N a bM a (y intercept) = 5.8158 b (slope) = 0.9726 b Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 39 The b-value represents the relative occurrence of large and small events in an event population. A high b-value means a higher proportion of low magnitude events to large ones and vice versa. Recent studies have also examined b values in microseismic datasets (e.g. Wessels et al. 2011) which unlike natural earthquake are often higher than 1.0. Many investigations, most of them using global or regional seismic catalogues, have studied the controls on b value. For example, b-value has been shown to correlate negatively with the normalized stress intensity (Hatton and Main, 1993). Schorlemmer et al. (2005) found that for global natural earthquakes, normal faulting and thrust events tend to have higher and lower b values, respectively. Higher b-values are expected for a complex fracture network and lower bvalues for planar fracture networks (Henderson et al. 1999). The relation between b values and different aspects of a failure regime is summarized in Table 5.1. Table 5.1. Expected influence of stress, fluids, source mechanism and fracture network complexity on seismic b value (Verdon, 2013). 5.1.1 Estimation of b-value In order to calculate b-values, a Matlab code created by James Verdon in 2012 (see Appendix B) has been used. The code determines seismic b-value from the gradient of log10NM against M. In order to estimate b-values accurately it is necessary to define a minimum magnitude (MMIN) above which all the events occurred in the volume of study can be detected. The code loops over increasing cut-off magnitude until a minimum magnitude (MMIN) at which the observed distribution can be modelled by the equation (7) is found. Then, using a Kolmogorov-Smirnov test, null hypothesis are rejected at a 20% significance level. If a MMIN at which the approximation of the observed data give a significance level lower than 20%, is not found, the b-value is not estimated. Stress Fluids Mechanism Network complexity b increase Low stress Fluid playing a role Tensile/extensional Complex network b decrease High stress Fluid not important Reverse Planar fracture/fault Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 40 -4 -3.9 -3.8 -3.7 -3.6 -3.5 -3.4 -3.3 -3.2 100 101 102 Moment magnitude Nº events (logarithmic scale) b-value estimation of Stage 1 event population Once MMIN is determined, the b value is estimated using the maximum likelihood method described by Aki (1965), where (5.2) and ⟨M⟩ is the mean of the magnitude distribution. To determine the error limits of b, MMIN is kept constant, and the values of b at which the Kolmogorov-Smirnov test is rejected at the 5% are used. The code is applied to the event population of every stage separately. Unfortunately, the event populations of the stages 2 and 4 are too small to estimate reliable b-values. Figures 5.2 and 5.3 show the results of the b-value estimation for stage 1 and 3 event populations respectively. The b-value of the stage 1 is higher than b-value of the stage 3, indicating that in stage 1 there were a higher proportion of lower magnitude events to larger ones than in stage 3. However, we should consider the estimation of the b-value for the stage 1 less reliable because of the difference between the maximum and minimum b-value limits. Figure 5.2. Approximation of the b value for the Stage 1 microseismic event population. The green and blue traces shows the magnitude distribution observed, with a moment magnitude lower and higher than the chosen MMIN, respectively. The red line shows magnitude distribution modelled by the equation (7) with the approximated b value. Finally, the dashed red lines are the error limits of the b value. a (y intercept) = -13.9398 best estimation of b (slope) = 4.2316 b bMIN = 2.6416 bMAX = 5.0316 MMIN = -3.6849 10 log MIN e bMM  Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 41 -5 -4.8 -4.6 -4.4 -4.2 -4 -3.8 -3.6 -3.4 -3.2 -3 100 101 102 Moment magnitude Nº events (logarithmic scale) b-value estimation of Stage 3 event population Figure 5.3. Approximation of the b value for the Stage 2 microseismic event population. See Figure 5.3 for graph explanation. Once b-values for both event populations are estimated, it could be interesting to see how bvalue varies spatially. Therefore, 2 by 2 meters x and y grids are applied to the different stage stimulation areas. A radius value is needed to define the microseismic event population of every point of the grid in order to estimate the correspondent b-value. The best results are obtained with a radius of 75 meters. These results are used to generate the b-values contour map (Figure 5.4). MMIN = -4.2646 a (y intercept) = -6.8884 best estimation of b (slope) = 2.0419 b bMAX = 2.2919 bMIN = 1.4619 Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 42 -50 0 50 100 150 200 250 Distance relative to the Observer well [ W (-) / E (+) ] (m) b values contour map (Magnitude-frequency distribution) -350 -300 -250 -200 -150 -100 -50 0 50 100 150 Distance relative to the Observer well [ S (-) / N (+) ] (m) 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 5 5.2 5.4 5.6 Figure 5.4. b-values contour map of the stages 1 and 3 with the location of the events, the observer wellhead and the geophones. Geophones Stimulation zone 1 (Stage 1) Stimulation zone 2 (Stage 3) Observer wellhead Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 43 5.2 Spatial distribution The spatial distribution for an event population can be described by the two-point correlation dimension (DC). The correlation integral, C(r), necessary to calculate DC describes the number of event pairs separated by a distance less than r: (5.3) where NP is the number of event pairs (NP=NE(NE-1)/2), NE is the total number of events and NP(R<r) is the number of event pairs separated by a distance less than r. If the events are distributed in a fractal manner, the two-point correlation dimension (DC) is related to the correlation integral by (5.4) therefore DC can be determined from the gradient of log10(C(r)) against log10(r). Generally, an event population can take three different forms: events following a linear feature, events delineating a planar feature and events distributed randomly through a 3D volume. Spatial distribution information that can be extracted from the two point correlation dimension DC has been demonstrated by Verdon et al. (2012) with a simple synthetic example (Figure 5.5). A DC ≈1 is expected for events distributed linearly, a DC ≈2 for events distributed in a plane and a DC ≈3 for events distributed randomly through a volume. Figure 5.5. Synthetic distributions for DC calculation: events distributed along a linear feature, events distributed on a plane, and events distributed randomly throughout a volume (Verdon, et al, 2013). 5.2.1 Calculation of DC A Matlab code also created by James Verdon in 2012 (see Appendix C) has been used to calculate the two-point correlation dimension (DC). The code calculates DC from the gradient of log10(C(r)) against log10(r). If the value of r increases beyond the extent of the event population, () C D C r r 2 ( ) () P P N R r Cr N   Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 44 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 Log10(r) Log10(C(r)) Dc calculation of Stage 1 event population the correlation integral C(r) will not have a significant increase, which implies a decrease of the gradient of log10C(r)/ log10(r). Therefore, similarly to MMIN in the calculation of b, a maximum cut-off for r, rMAX, must be defined. The code loops over decreasing cut-off r until a maximum r (rMAX) at which the observations fit a straight line is found. A residual misfit is computed to quantify the fitting of the model to the observations: (5.5) where Ciobserved is the observed C(r) at a given r, and Cimodelled is the corresponding modelled value. The rMAX is chosen as the maximum value of r at which MF ≤ 5%. To determine the error limits of DC, rMAX is kept constant and the values of DC at which the residual MF is 5% are used. The code has been applied to the event population of every stage separately. Figures 5.6, 5.7, 5.8 and 5.9 show the results of the DC value estimation for stages 1, 2, 3 and 4 event populations respectively. As the stage 4 only has 18 events, the DC value obtained might be distorted. The DC value for the stage 1 is higher than the DC values for the stages 2 and 3, with values of 2.5585, 2.3355 and 2.2409 respectively. Figure 5.6. Approximation of the DC value for the stage 1 microseismic event population. The green and blue lines represent log10C(r) against log10(r) for r values higher and lower than the rMAX respectively. The red line represents the C(r) values calculated with the approximated DC. Finally, the dashed red lines are the error limits of the DC value. Log10( rMAX = 63.8822 ) = 1.8054 best estimation of Dc (slope) = 2.5585 Dc DcMAX = 1.6285 DcMIN = 3.1485 mod 100 100 MAX i robserved eled observed i i i ri MF C C C        Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 45 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 Log10(r) Log10(C(r)) Dc calculation of Stage 2 event population 0 0.5 1 1.5 2 2.5 -4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 Log10(r) Log10(C(r)) Dc calculation of Stage 3 event population Figure 5.7. Approximation of the DC value for the stage 2 microseismic event population. See Figure 5.6 for graph explanation. Figure 5.8. Approximation of the DC value for the stage 3 microseismic event population. See Figure 5.6 for graph explanation. Log10( rMAX = 65.4375) = 1.8158 best estimation of Dc (slope) = 2.3355 Dc DcMAX = 1.4955 DcMIN = 2.8855 Log10( rMAX = 20.7357 ) = 1.3167 best estimation of Dc (slope) = 2.2409 Dc DcMAX = 1.4109 DcMIN = 2.7209 Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 52 Polarization of the fast wave (ψ) is expected to be parallel to fracture strike for vertical waves and parallel to sedimentary layering for horizontal waves (Figure 6.2). Therefore, shear-wave splitting measurements can be used to identify different anisotropy sources present in a rock by making some simplifying assumptions about their orientation and symmetry. As a sedimentary hydrocarbon reservoir is studied in this project, the anisotropy of the sedimentary layering is considered horizontal or subhorizontal such that SWS is maximum for horizontally propagating waves, whereas fracture sets are assumed to be vertical or almost vertical such that splitting is maximum for vertical waves. Additionally, the amount of splitting is proportional to fracture density which is a measure of how fractured the rock is. Figure 6.2. Schematical diagram showing shear-wave splitting for vertically and horizontally propagating waves in a medium with vertical fractures (top) and horizontal bedding (bottom) (Wuestefeld et al. 2010). Since the propagation of the shear-waves is not always vertical or horizontal, a more complex model is needed to predict the SWS results for various propagation directions. Upper hemisphere projection is often used to show SWS for different propagation directions obtaining a geographical visualization of the results (Figure 6.3). By measuring SWS along ray paths at various directions has been found that variations in fast shear-wave polarizations and delay times reflect the symmetry and the strength of anisotropy (e.g. Wuestefeld et al. 2010) (Figure 6.4). For example, in the Figure 6.4(a) horizontal bedding will produce anisotropy with hexagonal symmetry with a vertical axis of symmetry (vertical transverse isotropy, VTI) which Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 53 means that maximum anisotropy is obtained for horizontally propagating waves whereas no splitting is obtained for vertical propagating waves. In the Figure 6.4(b) almost vertical or vertical fracture sets produce horizontal transverse isotropy (HTI) which means that no splitting is obtained for horizontal waves propagating perpendicular to the orientation of the fractures whereas maximum anisotropy is obtained for waves, with any inclination, propagating parallel to the orientation of the fractures. Finally, the combination of these two anisotropy sources produces a more complex model that can be seen in the Figure 6.4(c). Figure 6.3. Schematic diagram showing how fast wave polarization (ψ) is deteermined in an upper hemisphere projection. Figure 6.4. Synthetic upper hemisphere plots showing SWS magnitude (colour contours), δVS (tick lengths), fast wave polarization, ψ (black tick orientations) and azimuth and inclination of the ray path for: (a) VTI anisotropy due to horizontal layering/fabric; (b) HTI anisotropy due to aligned vertical fractures and (c) orthorhombic anisotropy due to vertical fractures in a horizontally layered medium (Baird et al. 2013) Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 54 0 1 2 3 Sen.01 Sen.02 Sen.03 Sen.04 Sen.05 Sen.06 Sen.07 Sen.08 Sen.09 Sen.10 Sen.11 Sen.12 S-Az = 123 S-Az = 125 S-Az = 100 S-Az = 134 S-Az = 128 S-Az = 135 S-Az = 116 S-Az = 111 S-Az = 146 S-Az = 85 0.15 0.2 0.25 0.3 0.35 Sen.01 Sen.02 Sen.03 Sen.04 Sen.05 Sen.06 Sen.07 Sen.08 Sen.09 Sen.10 Sen.11 Sen.12 Time (seconds) 6.1 Shear wave splitting measurements The location of the events is needed for shear-wave splitting measurements. Hence, the event data of the acquisition company is used to provide a measure of anisotropy within the reservoir. Stages 1 and 2 are located southeast from the receivers whereas stages 3 and 4 are located northeast from them. It is important to remind that the direction of maximum stress is approximately northeast-southwest. Based on the locations of the stages some assumptions are made about which anisotropy source could be reflected in the SWS measurements of every stage. The ray path between the events and the receiver from stages 1 and 2 should be subhorizontal and approximately perpendicular to the maximum direction of stress. Therefore, we expect to obtain an estimation of the anisotropy, due to horizontal layering/fabric, from the SWS measurements from stages 1 and 2. The ray path between the events and the receiver from stages 3 and 4 should be steeply propagating and approximately in the orientation of maximum direction of stress. So an estimation of the anisotropy, dominated by aligned vertical fractures, is expected to be obtained from the SWS measurements from stages 3 and 4. Fractures are expected to be oriented in the direction of maximum stress, so we expect SWS measurements oriented roughly in the same direction. The source-receiver paths do not travel outside the reservoir rocks, therefore we assume that there is little spatial variation anisotropy. However, the splitting measurements will vary depending on the propagation direction of the shear wave. All the S-waves arrivals of the different stages are picked manually (Figure 6.5). Figure 6.5. Traces and S-waves picks of an event file (left graph); particle motion plot and direction of propagation (perpendicular to the particle motion orientation) for the S-waves picked (right graph). Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 55 Then, the data is analysed using the automated splitting approach of Wuestefeld et al. (2010). It rotates the traces into the ray frame (SH, SV, P), where P is the direction of propagation which should contain most P wave energy, SH is the horizontal component perpendicular to P and SV is the component perpendicular to P and SH. S wave energy should be partitioned between the components SH and SV. The method finds the splitting within the SH-SV plane. Splitting is estimated by windowing the S-arrival and doing a grid search through splitting parameters (ψ and δt), which are used to remove the splitting. The best estimate should linearize the particle motion in the SH-SV plane. This method is easily applied to large datasets and provides a quality control that consists on defining a quality index which varies from -1.0 for null measurements, to 0.0 for poor and +1.0 for good measurements. The quality index is based on differences between SWS estimations using two different techniques (eigenvalue and crosscorrelation methods). A diagnostic plot showing all this process is created for every splitting measurement. Diagnostic plot for good measurement from stage 1 can be seen in the Figure 6.6. We define good measurements as having a quality index higher than 0.8, a signal-noise ratio higher than 4, a time lag less than 4.5 ms with an error less than 0.5 ms and an error in fast polarization of less than 15º. Figure 6.6 Diagnostic plot for a good splitting measurement in the stage 1. The top left panel shows SH (red), SV (blue) and ray (black) seismogram components and S window (yellow). The top right panel shows radial (green) and transverse (magenta) component before (top two traces) and after (bottom) splitting correction. The lower left panel shows in the two top graphs the SH and SV waves before (left) and after (right) the splitting correction; in the two bottom graphs the particle motion in SH-SV coordinates before (red) and after (blue) correction (left) and the quality index for all the measurements (right). Finally, the lower right panel shows the error surfaces of the eigenvalue (left) and crosscorrelation (lower right) methods. These represent the error surfaces of the respective best measurements. The location of the minimums is indicated by thin black lines. All the fast direction measurements and their respective time-lag are plotted in the higher right graph. Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 56 6.2 Interpretation of the results Comparing the good shear-wave splitting measurements projected in upper hemispheres with the synthetic upper hemisphere plots in the Figure 6.4, it is possible to see the observed anisotropy matches our expectations. As we predicted, the shear-wave splitting measurements for stages 1 and 2 (Figure 6.7) are consistent with vertical transverse isotropy (VTI) as would be expected from horizontal bedding (Figure 6.4a). That can be seen clearly in the results of the geophone 5, where all the polarizations are subhorizontal. The anisotropy is generally low with an average of 4% which is reasonable for a horizontal layering/fabric. The Figure 6.8 shows upper hemisphere plots of the resulting splitting measurements for stages 3 and 4, which are consistent with horizontal transverse isotropy (HTI) as would be expected from vertical fracture sets (Figure 6.4b). Therefore, as expected the splitting measurements in stages 3 and 4 are likely related to aligned vertical fractures. Additionally, it can be seen that the orientation of almost all the fast polarized waves is similar to the maximum direction of stress orientation which means that the fractures are roughly oriented in the direction of maximum stress, as we predicted. In this case, the anisotropy estimated is higher with an average of 10%. Splitting measurements provide enough information to visualize the anisotropy symmetry of the medium, consisting of horizontal bedding and vertical fractures oriented in the direction of maximum stress. However, it is important to note that the medium anisotropy has been simplified a lot. If a more detailed estimation of the anisotropy is required, more data would be necessary. Unfortunately, insufficient high quality SWS measurements have been obtained to adequately constrain an inversion for detailed fracture and rock fabric properties in the rock surrounding the stimulated volume. However, inversion techniques have been used to infer quantitative fracture properties from larger datasets (e.g. Verdon and Wuestefeld 2012; Baird et al. 2013). Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 57 Figure 6.7. Upper hemisphere plots of the resulting splitting measurements for stages 1 and 2. The top hemisphere plot shows all the splitting measurements whereas the bottom ones show the splitting measurements separated by geophones. The plots show the azimuth and inclination of the S-wave arrival, the geographical orientation of the fast polarized wave (ticks) and the anisotropy (colour) calculated with the δt. Geophone 5 Geophone 3 Geophone 8 Geophone 2 Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 58 Figure 6.8. Upper hemisphere plots of the resulting splitting measurements for stages 3 and 4. See explanation of the plots in the Figure 6.7. Geophone 8 Geophone 5 Geophone 9 Geophone 4 Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 59 7. Conclusions and future research 7.1 Conclusions A new automated microseismic event detection method was developed and applied to a continuous microseismic dataset obtaining 1098 event detections. Cross-correlation was proved as a good technique to improve the event detection limits of the STA/LTA method, reducing number of false detections and increasing the confidence of the picks. Having the results of the acquisition company was essential in order to assess how effective the technique was. The data gaps identified in the event data provided by the acquisition company has been filled by the new method. Additionally, while the events detected by the acquisition company occurred during the second and third days of monitoring exclusively, almost half the events detected by the new method occurred during the first day of monitoring, providing more information for an assessment of a future project. Therefore, the new automated microseismic event detection method has produced useful results and presents an alternative to other automated methods developed previously. Since we did not have enough time to estimate the location and magnitude of the events detected by the new method, further analysis was done with the event data provided by the acquisition company. Magnitude and spatial distribution maps, composed by b and DC values respectively, have provided a good insight of the mode of failure and the complexity of the fracture network induced. The 4 different stimulation stages were studied separately. Stage 1 presents a complex fracture network, events occurring throughout a volume whereas stages 2, 3 and 4 presents a simpler fracture network, events occurring along an approximately planar feature. Various possibilities were contemplated to explain this discrepancy, such as different injection conditions, or the presence of different natural fracture networks prior to the stimulation. Unfortunately, detailed well logs and injection rates were not provided by the acquisition company. The possible failure of the packer between stages 1 and 2 could be a differential factor as well. Additionally, the errors associated with the estimation of locations and magnitudes of the event data provided by the acquisition company are not known. As stage 1 is farther away from the receivers we could assume that the location and magnitude values of these events are less reliable and hence, the more complex fracture network deduced for stage 1 may be partially a consequence of poor resolved estimations of the locations and magnitudes of the events. However, as the results for stage 2 are similar to those obtained for stage 3 and are as far away from the receivers as the events for stage 1, it is unlikely to think that only the Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 60 locations and magnitudes of the events for stage 1 are poorly resolved. It is clear that more information is needed to explain with more confidence the different results obtained for stage 1. Finally, shear-wave splitting was used to estimate the anisotropy of the medium. SWS measurements for stages 1 and 2 were likely related to the horizontal bedding showing an average anisotropy of 4% and mostly subhorizontal polarizations, conversely the SWS measurements for stages 3 and 4 were likely related to the vertical fractures, oriented in the direction of maximum stress, showing an average anisotropy of 10%. The fractures are the dominant anisotropy of the medium because they induce a larger effect in the shear-wave splitting. Unfortunately, it was not possible to constrain a detailed fracture inversion because of insufficient high quality SWS measurements. However, we consider that SWS measurements provided important information about the symmetry anisotropy and the orientation of the fractures. 7.2 Future research Future research could be done from the work of this project. It would be really interesting to develop an approach to estimate the locations and magnitudes of the microseismic events detected by the new method. With this information it would be easier to identify what activity produced the events detected on the first day of the monitoring. Then, it could be possible to do further analysis with the new event data. For example, a study of the focal mechanism of the events will provide information about the type of fault in the different stages, complementing the information about the fracture networks obtained in the present project. Actually, the long term objective is to create a common workflow for microseismic datasets acquired from stimulations of unconventional reservoirs, to obtain at the end, as much information as possible about the natural and induced fracture network. For example, this workflow could consist of applying the automated method developed in this project to detect the events, using an approach to estimate the locations and magnitudes of the events and doing further analysis such as magnitude and spatial distributions, focal mechanism of the events and shear-wave splitting. Development of an automated microseismic event detection method and analysis of a microseismic dataset from hydraulic fracture stimulation to characterise natural and induced fracture networks within a tight reservoir 61 8. References Aki, K., 1965. Maximum likelihood estimate of b in the formula log N = a – bM and its Confidence Limits. Bulletin of the Earthquake Research Institute, University of Tokyo, 43, p. 237-239. Akram, J. & Eaton, D., 2012. Adaptive microseismic event detection and automatic time picking. University of Calgary, Canada. GeoConvention 2012: Vision. Allen, R., 1978. Automatic earthquake recognition and timing from single traces. Bulletin of the Seismological Society of America, 68(5), p. 1521-1532. Baird. A.F.; Kendall, J.-M.; Verdon, J.P.; Wüestefeld, A.; Noble, T.E.; Li, Y.; Dutko, M. & Fisher, Q.J., 2013. Monitoring increases in fracture connectivity during hydraulic stimulations from temporal variations in shear wave splitting polarization. Geophysical Journal International, 195, p. 1120-1131. Bakken Decision Support System, Energy & Environmental Research Center (http://www.undeerc.org/Bakken/bakkenformation.aspx, last access 02 05 2014). Chambers, K.; Kendall, J.-M.; Brandsberg-Dahl, S. & Rueda, J., 2010. Testing the ability of surface arrays to monitor microseismic activity. Geophysical Prospecting, 58(5), p. 821-830. Forghani-Arani, F.; Behura, J.; Haines, S.S. & Batzle, M., 2012. An automated cross-correlation based event detection technique and its application to a surface passive dataset. Geophysical Prospecting, 61, p. 778-787. Garnero, E.J., Images relating to seismic anisotropy in Earth’s mantle (http://garnero.asu.edu/research_images/images_anisotropy.html, last access 29 04 2014) Hatton, C.G.; Main, I.G. & Meredith, P.G., 1993. A comparison of seismic and structural measurements of scaling exponents during tensile subcritical crack growth. Journal of Structural Geology, 15(12), p. 1485-1495.