Abstract
Despite the current easing in demand for increased oil production linked to the global downturn in crude prices, energy demand continuously increases and the long-term demand will require maximizing the productivity of reservoirs and a search into the exploitation of new resources in increasingly challenging environments. In this study, we present the results from the monitoring of the very first multistage stimulation experiment at a shale gas reservoir in Saudi Arabia, presenting an analysis of the microseismicity induced during the treatment. Our aim was to analyse microseismic events to better understand fracture growth and the role of pre-existing fractures in these reservoirs. Microseismic (MS) event monitoring is used to track the creation of fractures during and after the stimulation, and therefore to evaluate the effect of the reservoir stimulation. The monitoring includes a downhole array of 12 3C-sensors that were deployed in a vertical well with a 30.5 m level spacing. A total of 415 MS events were located and analysed, with the results outlining induced fractures extending consistently with an average azimuth of N335° E, normal to the horizontal section of the treatment well. This implies that there are no changes in the local stress direction along the treatment well either in situ or induced along the treatment. There are significant changes in total length and aspect ratio (length/width) of the fractures induced in the different stages. These variations could be attributed to in situ fracturing, local rock heterogeneity or the influence of the treatment parameters. In general, early and late stages of stimulation show the longest fracture networks, with events induced further away from the initiation point. We found no immediate relationship between treatment parameters (peak pressure and pumping rates) and fracture extension. Sensitivity analysis using Monte Carlo simulation methods shows a higher location uncertainty for events located at the early stages, thus limiting the interpretation from monitored seismicity in the early stages. An analysis of magnitude distribution with distance shows a decrease in sensitivity of one degree of magnitude for every 375 m, and a maximum viewing distance of approximately 700 m for the current set-up. The low number of located events does not provide a complete enough dataset for a robust analysis of changes in b-value (slope in linear part of magnitude distribution) during the treatment: however, magnitude distributions, corrected for array sensitivity, provide a useful variable for the validation of geomechanical models currently being developed for the reservoir.
Globally, shale gas resources represent nearly two-thirds of the world's hydrocarbon reserves. The United States, for instance, produces shale gas in no less than 15 States and production has remarkably increased during the last 10 years (Duggan-Haas et al. 2013). According to the United States Energy Information Administration, almost 50% of gas production in the USA in the next 20 years will be likely to come from shale gas (EIA 2013).
Since shale gas is trapped within a low permeable layer, it is challenging to produce economically and, in most cases, hydraulic fracturing (HF) technology is required to stimulate production. The HF technology involves the injection of high-pressure fluid, proppant and small amount of chemicals to create fracture networks enabling the trapped gas to flow towards the production well, thereby increasing conductivity and gas flow (Green et al. 2007). Thus, during stimulation, it is necessary to monitor the fracture growth and direction within the shale gas layer. Relatively recently, various technologies have been used to image the hydraulic fracture growth. Among them is microseismic monitoring, which has become an important technology for imaging the evolution of treated reservoirs, providing near-real-time feedback on the propagation of the induced fractures (e.g. De Meersman et al. 2009; Verdon et al. 2010). In fact, during the last 10 years, microseismic monitoring technology has made significant contributions in imaging the geometry of fractures associated with hydraulic stimulation (Maxwell et al. 2010).
One of the most important steps in microseismic monitoring is to detect, locate and, more importantly, to analyse microseismic events. Analysing microseismic events is also useful for interpreting fracture location, orientation and extent in all directions of the stimulated fracture network. This enables the effect of the HF treatment to be tracked in near-real-time and the extent of stimulated volume to be estimated (e.g. Warpinski 2009; Maxwell et al. 2010). Consequently, we have undertaken this study to analyse microseismic events to better understand fracture growth in this new scenario. This study was carried out as part of a collaborative effort between Saudi Aramco and the King Fahd University of Petroleum and Minerals for the development of reservoir stimulation in the Arabian Peninsula region. The collaborative effort is centred on developing a methodology for generating dynamic, high-resolution seismic and geomechanical models of shale reservoirs before, during and after stimulation, and interpreting the models in terms of fracture susceptibility and fracture dynamics. We use the general seismic cloud geometry and magnitude distribution as the observational validation tool for geomechanical models, being more robust against uncertainties than the location and characteristics of individual events. Figure 1 shows an overview of the microseismic analysis workflow, which was followed in this study to validate predictive numerical models for the stimulation of a shale gas reservoir. This study is the first in the region, in the NW of Saudi Arabia, to analyse the first multistage hydraulic fracture stimulation programme carried out to target the lower Silurian Qusaiba shales, and therefore the response observed in this study is key for future developments of shale gas resources from this formation.
Schematic of the workflow followed in this study for the analysis of induced microseismic activity to validate predictive numerical models for the stimulation of a shale gas reservoir.
Geological setting
The Qusaiba is a Lower Silurian clastic assemblage and forms part of the Qalibah Formation (Jones & Stump 1999). On a regional basis, the Lower Silurian clastics represent a progradational sequence that was deposited on a broad continental shelf of Gondwana during the sea-level rise following the deglaciation of the Gondwana polar ice cap (Bender 1975). The Qusaiba consists primarily of well-bedded to laminated alternations of mainly grey to black micaceous shale and siltstone, with thin sandstone intercalations and micaceous partings. The lower few tens of metres of the Qusaiba in the NW of Saudi Arabia is informally identified as the ‘hot shale’, because it exhibits a diagnostic sharp increase in gamma radiation typical of many organic-rich facies. Occasionally, above the ‘hot shale’, are very thin black shales that show similar log characteristics (Mahmoud et al. 1992). The shale facies were generally deposited under dysoxic to anoxic conditions in an intra-shelf basin located north of the Arabian Arch (Fox & Ahlbrandt 2002). The Qusaiba disconformably overlies the Sarah Formation, and the contact generates a regionally mappable, high-amplitude seismic reflection. Silurian ‘hot shale’ source rocks began generating hydrocarbons possibly during the late Palaeozoic in the region of northern Arabia into Iraq, and during the Triassic in the region of the present-day Arabian Gulf (Bishop 1995).
Present-day maximum horizontal in situ stress in the area is dominated by Arabian plate tectonics, with a trend centred around the NW–SE direction. Borehole image logs and core data show fractures from micro- to macroscales, including extension fractures and compressional shear fractures (Ameen 2014). Ameen (2014) also noted that forward modelling suggests that the natural fracture networks are not critically stressed under reservoir conditions, but when subjected to hydrofracture stimulation they and bedding discontinuities could form the seeds for a complex hydrofracture network.
Dataset for the study
Our dataset includes microseismic events recorded through a comprehensive monitoring programme initiated by Saudi Aramco for HF operations. This monitoring programme was the very first multistage stimulation experiment for unconventional shale gas reservoirs in the Arabian Peninsula. The monitoring programme consists of a downhole array, with sensors cemented in shallow wells and a number of surface arrays. In this study, we only processed data from the downhole array. Surface array data were characterized by low signal to noise ratios (SNRs) such that only the very largest subsurface hydrofracture events could be detected, such as the perforation shots and events created during the initial fracture generation. Hence, the results from analysis of the surface monitoring data were not useful in determining any details about the more numerous, weaker events and were not included in this study. The downhole array consists of 12 three-component sensors deployed along a neighbouring vertical well with an aperture of 335 m centred on the target level. The HF treatment was carried out along the horizontal section of a deviated well through 12 stages, performed over a 10 day period, using the plug-and-perf method.
Data processing
We first removed a constant voltage offset (i.e. DC offset) that appears on the waveforms due to the acquisition electronics. The removal is performed by calculating the mean amplitude in a window of waveform data. This window is taken on the constant signal preceding the waveform arrivals associated with the event and is taken with sufficient length to average out any noise.
We then applied a bandpass filter between 20 and 600 Hz to enhance the SNR of the traces.
The next step was the processing of the continuous data stream to harvest and locate microseismic events using a block velocity model based on a VSP survey carried out at the monitoring well (Fig. 2). Owing to the low resolution of the survey, only seven layers were defined, with two of them at the level of the monitoring tools.
Block velocity model (dashed lines) used in the processing of the microseismic events in this study built from vertical seismic profile (VSP) surveys (solid lines).
Event locations were performed using a grid-search optimization approach, performing an exhaustive search over the complete reservoir for the minimum in an objective function that includes differences in travel times for P- and S-wave arrivals, and dip and azimuth for source vectors, using an L2 norm for the misfit calculation. In order to optimize the processing time and resolution, a collapsing approach was used, with an initial search over a coarse grid of 30.5 m spacing reducing to a final resolution of 10 m. Forward modelling was performed using a two-point ray-tracing algorithm proposed by Tian & Chuen (2005). The algorithm includes the calculation of source vectors and travel time for direct and reflected body waves and any potential head-waves generated on high-contrast interfaces. In this case, no head-waves had been observed or were estimated as first arrivals from the model.
We scanned the full data stream to extract discrete waveform sections corresponding to potential microseismic events that included both P- and S-wave arrivals. The search was performed using a double window across the square root of the sum of squares of the full microseismic record, defining a triggering function F for each waveform data point i for each trace as:
(1)where Aj is the amplitude, NSW is the length of the front search window or signal window in data points, and NNW is the length of the back search window or noise window in data points. On average, values of 400 and 200 ms for the back and front windows, respectively, were used for event triggering, with wider windows (500 and 250 ms) for early more distant stages. Peaks in the function indicate a sudden increase in amplitude relative to the data behind it: this difference is used as a threshold criterion for defining a trigger. An event is expected to have coherent arrivals recorded by multiple channels within a time window: therefore, further criteria were set for a minimum number of tools to exceed the threshold in the trigger function within a set time window. This time window is dependent on the tool spacing, while the section of waveform that is extracted corresponding to a potential event (trigger) depends on the distance between the monitoring array and the well stimulation locations. Owing to the extent of the stimulated area, we used two different waveform lengths: 1 s for early stages, 1 – 5; and 0.5 s for late stages, 6 – 12, positioned closer to the monitoring well.
For each located event, an estimate of the released energy is provided by a location magnitude (ML) calculated using the RMS of the velocity waveform amplitude for a window around the P- and S-wave arrivals for the instruments and phases used in the location of the event:
(2)where WRMSi is the mean RMS velocity amplitude for the three traces in instrument i, and di is the length of the ray path between the event and instrument i. We note that moment-magnitude is the most common definition used: however, using local magnitude as an event size estimate is sufficient for the purpose of this study.
Following the event scanning, we observed a large number of ‘tube waves’ that showed coherent large-amplitude arrivals in all monitoring tools, but with an unrealistically low apparent propagation velocity. Examples of the harvested waveforms are shown in Figure 3. These waves have been observed in other hydraulic treatments, and are interpreted as having been generated from the perturbation of fluid in the wellbore, travelling up and down the wellbore, along the interface between the fluid in the wellbore and the wall.
Example waveforms for triggers harvested from the continuous microseismic record during the stimulation of the shale gas reservoir including: (a) an event recorded during the stimulation of Stage 6 corresponding to a ‘tube wave’ induced at the wellbore interface; and (b) waveforms for a located microseismic event recorded during the stimulation of Stage 10.
In general, the SNR for the microseismic (MS) events analysed in this experiment was relatively low, with a median value of 1.66, which required manual checking of identified phases in most cases. The number of tools in the array allowed for a robust location in most cases, after discarding hard to identify outlier events.
Results
We successfully located a total of 415 MS events during the entire reservoir stimulation (Fig. 4). Figure 5 shows the locations for all the located events for the 12 stimulation stages, describing fractures extending consistently at an average azimuth of N335° E, normal to the horizontal section of the treatment well.
Histogram of the number of detected and located microseismic events along the 12 stages of the stimulation analysed in this study.
Location of the 415 events located during the hydraulic stimulation colour coded by time (dark colours correspond to early events, light colours to late events) and size scaled to relative magnitude in (a) the map view and (b) the side view looking north. Grid spacing is 50 m.
There were significant changes in total length and aspect ratio (length/width) of the fractures induced in the different stages that could be attributed to in situ fracturing (Ameen 2014), local rock heterogeneity or the influence of the treatment parameters. In general, early and late stages (2, 4, 8 and 10–12) show the longest fracture networks, with events triggered further away from the initiation point.
Since we used a low-resolution 1D velocity model for the location of the harvested events, it is necessary for us to analyse the impact of the poor determination of the velocity model on event locations. A typical approach for assessing the uncertainties in event location due to an unknown variability in multiple variables is running a Monte Carlo simulation (e.g. Billings et al. 1994). In this case, to analyse the impact of the uncertainties in the velocity model on the observed location results, the values of the P- and S-wave velocities were altered randomly following a Gaussian distribution for a fixed number of layers for the location inversion, repeating the process at two locations representing early and late stages. The resulting scatter in locations is plotted in Figure 6, showing similar vertical scatter for both test locations, while there is a significantly higher sensitivity in the horizontal uncertainty for events located at the early stages, with a scatter doubling in magnitude that of late stages. This analysis does not take into account other important sources of location uncertainty, such as errors in phase picking or uncertainties in tool orientation: however, with the monitoring geometry considered here, uncertainties in velocity models are the main source of location error (Jones et al. 2014).
Contour density plots showing the results from the Monte Carlo simulation to analyse the sensitivity of event location against uncertainties in the velocity model. Location inversion was performed 200 times at each theoretical position using a velocity model in the inversion with a 5 and 10% deviation from the theoretical model. The plots show the resulting event locations from each of the simulations; the results for simulations with 5% deviation are shown in (a) map view and (b) side view looking north, and the results for simulations with 10% deviation are shown in (c) map view and (d) side view looking north. The stars (yellow in the online version) mark the true locations of theoretical events located at an early stage (eastern/right-hand side of the plots) and at a late stage (western side of the plots). Dark grey events correspond to the resulting locations for the simulations at the early position, and light grey events correspond to the locations for the simulations at the late-stage position. Grid spacing is 30 m.
Located events are used as a first approximation to estimate the volume of the stimulated fracture zone, which is one of the main variables to be considered in the correlation and validation of geomechanical models. We estimated the stimulated volume in each stage assuming a Gaussian distribution of the events around each centroid and fitting the ellipsoids defined by the eigenvalues of the covariance matrix encompassing 90% of the events located during each of the stages (Fig. 7). An overlap was observed between the events observed during some of the treatment stages. This overlap was observed between stages 4 and 5, and, more significantly, between stages 6 and 9.
Ellipsoids fitted to the located microseismic events induced in each treatment stage to determine the stimulated volume and fracture network geometry for each stage.
No events were located for stages 1 and 3, where triggers corresponded mostly to tube waves or events with a very low SNR. Located events show a good convergence for the search for minima in the objective function based on travel time and source vector residual distribution for the complete, with median location residuals of 2.5 m. Note that this residual does not represent the uncertainty in event location as no potential errors (e.g. introduced by poor determination of the velocity model or errors in phase picking) are considered. The structure and geometry of the cloud of triggered events is more robust against location errors than individual events, and the fracture network geometry and relative changes in stimulated volume and aspect ratio between stages will be used as the primary information for correlation with numerical models.
Temporally, all located events show a good correlation between injection rates and triggered events, with activity initiating on the first step in pumping rate, increasing at each stage and in the case of late stages (stages 7–10), with activity extended for a period of 5–10 min after pumping was stopped, possibly associated with leak-off (Fig. 8).
Engineering treatment parameters compared to observed and located events for two example stages: Stage 5 (a) and Stage 7 (b). These stages showed different propagation distance despite similar treatment parameters being used across the stages.
The evolution of the ratio between the energy content of the S- and P-waves of the microseismic events can be used to gain some insight into the dominant mechanism of microfracturing, and was used as an estimate of the relative shear or tensile component in the rupture that originated the microseismic event (Fig. 9). A critical value of ES/EP < 10 (where ES and EP are the energy contents of the S- and P-wave phases, respectively) was taken as an indication of the existence of a significant tensile (volumetric) failure component (Gibowicz & Kijko 1994).
Evolution of the ratio ES/EP for all the located events.
An estimate of the sensitivity of the array and the completeness of the dataset for each stage is provided by the magnitude distribution with distance. As mentioned above, in this case a local magnitude scale (equation 2) is used to estimate the size of the event. Figure 10a shows that this monitoring configuration provided a maximum viewing distance of approximately 700 m, with a drop in sensitivity of one order of magnitude every 375 m. Figure 10b presents the typical magnitude distribution against the logarithm of the number of events for the calculation of b-values as the slope of the line fitting the decay in number of events for larger magnitudes (Fig. 9).
(a) Location magnitude (ML) v. distance: the dashed lines indicate the fitted sensitivity decay and the upper limit for the recorded magnitude. (b) Event magnitudes shown following the frequency–magnitude distribution developed by Gutenberg & Richter (1942) describing the number of seismic events (N) occurring above a certain magnitude (M) for all located MS events.
Discussion
The consistency in the direction of the propagation of MS events across the different stages indicates that there are no changes in the local stress direction along the treatment well either in situ or induced during the injection. The aspect ratio of the fractured volume in each stage presents significant changes across the stages, with fractures in mid stages showing a significant shorter propagation from the injection point. We find no immediate relationship between treatment parameters (peak pressure and pumping rates) and fracture extension that could explain these changes, as treatment parameters are similar for all stages (Fig. 8).
By only using microseismic information, the calculation of stimulated volume does not take into account any induced ‘smooth’ slip on joints or aseismic deformation, which must be considered when interpreting and correlating with the models. However, the relative changes in aspect ratio can be assumed to provide a valid correlation variable, providing a relative comparison of the effectiveness of the treatment at each stage.
An overlap is observed between the events located during stages 4 and 5, and, more significantly, between stages 6 and 9. This overlap indicates a degree of interaction between the induced fracture networks, the presence of in situ fractures or is possibly due to location uncertainty.
The slope of the fitted line to the larger monitored magnitudes in the event magnitude distribution shows a high b-value (2.3), typically associated with induced seismicity, with a dominance of small fractures indicating a relatively small role of larger in situ fractures in the monitored seismicity. The low number of events successfully processed and the limitations of the monitoring system do not allow for a robust interpretation of this parameter or the analysis of its temporal evolution: however, the relative magnitude distribution can be used in the correlation, with synthetic seismicity generated by numerical models.
Energy partition between the P- and S-wave arrivals could provide an insight into the dominant fracture mechanism. The temporal evolution indicates an increasingly presence of a tensile component: however, in this case, it is interpreted as an effect of the higher attenuation subject by the P-wave component of the waveform, which becomes increasingly dominant as events locate closer to the monitoring array. Further analysis is needed here.
Conclusions
We have successfully processed and analysed the microseismic events generated during the very first multistage stimulation study in the Arabian Peninsula with the objective of gaining a better understanding of the fracture growth process in this context.
Microseismic events provide an insight into the heterogeneity within the treated volume, changes in the effectiveness of the HF and changes in local stress. The SNR of the triggered events is low, and the majority of the events required manual picking of the P-wave arrival. Energy partitioning between P- and S-wave phases indicated a dominant shearing component in all located events. The energy content in the P-wave component increased as stimulation progresses. This indicates that the relative depletion in P-wave energy is an effect of the higher attenuation of this component, with higher dominant frequencies. The located microseismic events describe a logical pattern, with initiation points for each stage progressing towards the monitoring well as stimulation progress. The lineation described by the initiation points correspond with the direction of the treatment well, validating the accuracy of the processed locations.
We observed that microseismic activity increased towards the end of pumping, and continued for up to 10 min after the pump shut down in the mid and late stages. These shut-in events are likely to be triggered by the increase in pore pressure at the edge of the stimulated zone. We also observed an overlap between the events observed during stages 4 and 5, and, more significantly, between 6 and 9, indicating a degree of interaction between the induced fracture networks. We found that induced fractures extended consistently with an average azimuth of N335° E in all stimulation stages, indicating that there are no changes in the local stress direction. A possible indication of the reservoir in situ fracturing is provided by the marked changes in total length and width of fractures from stage to stage, implying the heterogeneity along the treated volume.
This analysis will serve as the basis for correlation with numerical models to determine what role different variables have on fracture propagation patterns. This will allow validation of the models for better prediction for the stimulation of shale gas reservoirs in this new geomechanical context.
Acknowledgements
The cooperation and invaluable assistance given to us by the Geophysics Technology Team at the EXPEC Advanced Research Center are gratefully acknowledged. We thank KFUPM for its continuous support of our research efforts. We also acknowledge the editor and anonymous reviewers for their valuable comments and suggestions that helped us to improve the final version of the manuscript.
Funding
This study was supported by Saudi Aramco under project number ES002357 at the King Fahd University of Petroleum and Minerals (KFUPM).
- © 2017 The Author(s)
References
Analysis of microseismic events during a multistage hydraulic stimulation experiment at a shale gas reservoir
