Received 14 April 2016; accepted 8 July 2016; published 11 July 2016
Hydraulic fracturing has become a common technology to enhance oil and gas production from low permeability rocks. Numerous microseismic events are generated during a fracturing process, and microseismic imaging has been used to monitor fracture propagation and development. Knowing how fractures grow helps engineers to determine optimal drilling directions and the number and spatial interval of fracturing stages  . Mapping microseismicity is also important in estimating reservoir volume and production  .
The application of microseismic images relies on the accuracy of event locations, which depend on the quality of data and velocity models. In  Eisner L., M. Thornton and J. Griffin found that even 5% change in velocity model could result in significant biases in the locations of microseismic events. Velocities from sonic logs are often different from horizontal velocities calibrated using perforation shots  , indicating a common presence of anisotropy in reservoir rocks. In addition, Grechka, V., and A. A. Duchkov concluded in  that hydraulic fracturing could cause remarkable changes in anisotropy. Therefore, ignoring anisotropy must cause errors in microseismic event locations.
Several studies have conducted in assessing the effect of anisotropy in microseismic locations using analytical forward modeling method for VTI (Vertical Transverse Isotropy) media   or homogeneous anisotropic models  . However, seismic anisotropy in fractured sedimentary rocks could be more complicated than VTI structure or homogeneous layering due to the superposition of intrinsic foliation, horizontal bedding, and aligned cracks. Anisotropy not only affects ray paths and travel times of seismic waves but also influences seismic amplitude and polarization   . Therefore, in this study we quantify the biases in microseismic event locations using synthetic waveform data from different configurations of anisotropy.
2. Synthetic Data
We use a pseudo-spectral numerical modeling method  to generate synthetic waveforms recorded at two downhole arrays in a three-layer velocity model. This is the first-time application of this method to microseismic settings. It generates reliable waveforms from 3D heterogeneous and anisotropic velocity models with realistic source parameters. Seven different anisotropic configurations are considered here. We use NonLinLoc, a nonlinear, probabilistic, global-search earthquake location algorithm   to determine the hypocenters from picked synthetic travel time data. The obtained event locations are compared with the true locations to assess the impacts from different anisotropy configurations.
Figure 1 shows the layout of the velocity model in the numerical simulations. It has a dimension of 3.5 × 3.5 × 2.0 km in x (north-south), y (east-west), and z (vertical) directions. The bottoms of the three layers are at the depths of 0.5 km, 1.2 km and 2.0 km, respectively. Two vertical monitoring arrays (A1 and A2) are set at (x = 0.994 km, y = 0.504 km) and (x = 2.996 km, y = 2.002 km). Each array contains 10 receivers at a 15 m interval. These receivers collect three-component seismograms from three sources (S1, S2, and S3) all in the middle layer (Figure 1). The acquisition geometry follows the description in  that a monitoring well is typically 160 - 550 m away from the fracture site. A double-couple source with a vertical fault plane striking in NS or EW is used for all three events. Seven velocity models, including isotropy, 5%, 10%, and 15% VTI and HTI (Horizontal Transverse Isotropy) media, are used in generating synthetic waveforms.
An example of vertical synthetic seismograms recorded by array 2 (A2) from event 3 (S3) is displayed in Figure 2 for a 15% VTI model. P wave can be identified easily despite the small amplitude. The waveforms also
Figure 1. Geometry of model, sources (red stars), and receivers (black triangles).
Figure 2. Vertical component synthetic seismograms recorded by array 2 (A2) from event 3 (S3) for 15% HTI case.
show clear shear wave splitting, which is expected for the given anisotropic velocity model and receiver-source configuration. In addition, the waveforms show reflected shear waves from the base of the middle layer and from the surface arrive in sequence with different move-out trends. The arrivals of P and fast S waves are picked from the synthetic seismograms and used in locating the events.
3. Location Method
The location method we used, NonLinLoc  , conducts non-linear, probabilistic, and global search for hypocenters in a 3D isotropic model. Following the inversion approach of Tarantola A. and Valette B. described in  , the program produces a misfit function, an estimate of the posterior probability density function (PDF) for the spatial hypocenter location (x, y, z). The location having the minimum misfit (or maximum likelihood) of the complete PDFs is selected as the “optimal” hypocenter. The origin time of the event is then calculated from the determined event location and observed arrival times.
We apply the NonLinLoc to all synthetic data from seven anisotropic models. Effective isotropic models are created for the corresponding anisotropic models in using the NonLinLoc with a grid interval of 50 m. The workflow of NonLinLoc method (Figure 3) includes four major steps: 1) build a 3D velocity grid from the given velocity description; 2) generate the grid for travel time calculations from the 3D velocity model; 3) predict travel times for given hypocenter locations; 4) find the optimal earthquake locations in the 3D model.
4.1. Results from Isotropic Media
Given the source and receiver locations (Figure 1), only the top two layers are useful in the calculation. The input velocities and densities in the numerical simulation are Vp = 2.4 km/s, Vs = 1.2 km/s, ρ = 2.0 g/cm3 for the top layer, and Vp = 2.8 km/s, Vs = 1.4 km/s, ρ = 2.2 g/cm3 for the second layer. The picked P and S arrivals at both arrays are used to determine the locations of 3 events simultaneously. The input velocities for the NonLinLoc inversion are the same as those used in numerical waveform modeling.
The location results are shown in Figure 4 and Table 1. The determined event locations (blue) align well with the true locations (red) in this case. The maximum error is 55 m in z direction for event S1 (top), 54 m in y direction for event S2 (middle), and 41 m in z direction for event S3 (bottom). The misfit of the determined and
Figure 3. Workflow of NonLinLoc method.
Figure 4. True (red symbols) and determined hypocenters (blue symbols) for the isotropic model. Star, circle, and diamond are for event S1, S2, and S3, respectively.
Table 1. Hypocenter locations (the unit is km).
true hypocenter locations are 58 m for S1, 79 m for S2, and 41 m for S3. The average location error for the isotropy case is 59 m, which is close to the grid interval used in the search. The error could also partially arise from the uncertainties in phase picks and partially from the difference of wave propagation methods used in the numerical modeling and event location.
4.2. Results from VTI Media
We generate synthetic waveforms for 3 VTI models with anisotropy strength of 5%, 10%, and 15%, respectively. The model can be described by 5 independent elastic coefficients, C11, C12, C13, C33, and C44. And C66 can be calculated from C11 and C12. The values of these elastic coefficients are based on laboratory data for sedimentary rocks   and are shown in Table 2 for the top two layers. The density is given as 2.0 g/cm3 for the top layer and 2.2 g/cm3 for the middle layer. The strength of anisotropy is calculated as the difference between the fast (horizontal) and slow (vertical) P wave velocity relative to the fast velocity,. Isotropic models that effectively represent the three anisotropic models are used in locating the events.
Figure 5 shows the determined and true hypocenters for all three events and the values are also listed in Table 1. There is no systematic error distribution in a preferred direction along the three axes due to the trade-off of all model parameters (x, y, z for 3 events). The average location error of 3 events varies from 158m for 5% VTI, to 237 m for 10%VTI and 258 m for 15% VTI case. All of them are much greater than that for the isotropic model (59 m). These large misfits are largely due to ignoring anisotropy in velocity models.
4.3. Results from HTI Media
We generate the HTI models by rotating the 3 VTI models. The 5 independent elastic coefficients are changed accordingly and recorded in Table 2. The anisotropy strength in the 3 HTI models is still 5%, 10% and 15%, respectively. Synthetic waveforms are produced from all 3 events for each model and the P and S arrivals are used to locate the events.
The results are shown in Figure 6 and in Table 1. The average location errors of 3 events are 156 m for 5% HTI, 244 m for 10% HTI and 265 m for 15% HTI, similar to those for the corresponding VTI cases. They are all much larger than the error of 59 m for the isotropic model, reflecting the influence of anisotropy.
Table 2. Elastic coefficients Cij for VTI and HTI media.
Figure 5. Same as in Figure 4 but for 5% (left), 10% (middle), and 15% (right) VTI models.
Figure 6. Same as in Figure 4 but for 5% (left), 10% (middle), and 15% (right) HTI models.
The above experiments show that event locations can be determined well with the errors as small as the grid interval using the NonLinLoc method for correct isotropic velocity models. Using the effective isotropic models instead of the true anisotropic models in the inversion, the variation of errors in x, y, and z directions does not show a systematic pattern with changing anisotropy, which is due to the trade-off among all model parameters. As expected, we observe that the biases in hypocentral locations increase with increasing anisotropy strength. As the anisotropy strength changes from 5% to 10% and 15%, the average misfit enlarges from 158 m, to 237 m and 258 m for the VTI models and from 156 m, to 244 m and 265 m for the HTI models. These errors are up to 3 - 5 times of the grid interval in the velocity model and are clearly caused by ignoring anisotropy in the inversion. We have tested the contribution to the errors from the initial event locations. The test shows that the error changes about 10 m for the 200 m difference in initial locations, which is negligible compared with the large errors due to ignoring anisotropy. All these results suggest that even ignoring 5% anisotropy in microseismic imaging can cause significant errors in hypocentral location, more than 3 times of the error for isotropic models.
The effective isotropic velocity used in event location is between the horizontal and vertical velocity for VTI or HTI structure. For example, in 15% VTI case, the effective isotropic P wave velocity in the top layer is 2690 m/s while the vertical and horizontal velocities are 2467 m/s and 2908 m/s. In practice, velocity model is developed from sonic logs that measure vertical velocity or calibrated from perforation shots that tend to be more influenced by horizontal velocity as in  and  . The discrepancy between the true and determined locations should be larger than the values in our results if either vertical or horizontal velocity is used in locating the events. Even using the effective isotropic velocities, biases in event locations are significant for even 5% anisotropy. As anisotropy can vary during a hydraulic fracturing process that constantly creates new fractures, it is critical to solving anisotropy from microseismic data and use the determined anisotropy model to locate events at each stage.
We conduct experiments in locating microseismic events from synthetic data that are generated using different VTI and HTI anisotropic models. Our results show that event locations can be determined with the errors as small as the grid interval for correct isotropic velocity models. When the effective isotropic models are used for data from anisotropic models, the location errors become 3 - 5 times of the grid interval. The biases in hypocentral locations increase with increasing anisotropy strength. As the anisotropy strength changes from 5% to 10% and 15%, the average misfit enlarges from 158 m, to 237 m and 258 m for the VTI models and from 156 m, to 244 m and 265 m for the HTI models. The results suggest that even ignoring 5% anisotropy in microseismic imaging can cause significant errors in hypocentral location, more than 3 times of the error for isotropic models.
We thank Mr. Jay Krishnan for offering technical support on Linux clusters computation issues.
 Maxwell, S. (2014) Microseismic Imaging of Hydraulic Fracturing: Improved Engineering of Unconventional Shale Reservoirs. SEG Distinguished Instructor Courses.
 Eisner, L., Thornton, M. and Griffin, J. (2011) Challenges for Microseismic Monitoring. SEG San Antonio 2011 Annual Meeting, San Antonio, 1519-1523.
 Warpinski, N. (2009) Microseismic Monitoring: Inside and Out. Journal of Petroleum Technology, 61, 80-85.
 Grechka, V. and Duchkov, A.A. (2011) Narrow-Angle Representations of the Phase and Group Velocities and Their Applications in Anisotropic Velocity-Model Building for Microseismic Monitoring. Geophysics, 76, WC127.
 Bayuk, I.O., Chesnokov, E. and Ammerman, M. (2009) Why Anisotropy Is Important for Location of Microearthquake Events in Shale? SEG Houston 2009 International Exposition and Annual Meeting, Houston, 1632-1636.
 Zhou, Z., Ngoc, N., Long, T., Borthwick, T., Zhao, K. and Teng, K.H. (2011) TTI/HTI Anisotropy Fro Fracture/Faults Image inside Granite Basement. EAGE Vienna 2011 Conference & Exhibition, Vienna, 15221.
 Chen Z., Cox, B. and Perkins, C. (2013) The Influence of HTI Anisotropy on Microseismic Event Location—A Case Study from a Tight Gas Field in British Columbia. SEG Houston 2013 Annual Meeting, Houston, 2003-2007.
 Hung, S. and Forsyth, D.W. (1998) Modeling Anisotropic Wave Propagation in Oceanic Inhomogeneous Structures Using the Parallel Multidomain Pseudo-Spectral Method. Geophysical Journal International, 133, 726-740.
 Lomax, A., Michelini, A. and Curtis, A. (2009) Earthquake Location, Direct, Global-Search Methods, in Encyclopedia of Complexity and System Science, Part 5, 2449-2473, Springer, New York.
 Lomax, A., Virieux, J., Volant, P. and Berge, C. (2000) Probabilistic Earthquake Location in 3D and Layered Models: Introduction of a Metropolis-Gibbs Method and Comparison with Linear Locations. In: Thurber, C.H. and Rabinowitz, N., Eds., Advances in Seismic Event Location, Kluwer, Amsterdam, 101-134.
 Thomsen, L. (1986) Weak Elastic Anisotropy. Geophysics, 51, 1954-1966.
 Wang Z. (2002) Seismic Anisotropy in Sedimentary Rocks. Part 2—Laboratory data. Geophysics, 67, 1423-1440.