OR/17/002 Work Package 1: Development of monitoring tools for the characterisation and measurement of induced pressure and geomechanical changes in the reservoir

Chadwick, R A. 2017. DiSECCS Final Summary Report. Work Packages 1 – 4 British Geological Survey Internal Report, OR/17/002.

In this work package, we integrated hydromechanical simulation results with rock physics models and full-waveform seismic modelling to assess time-lapse seismic attributes for dynamic reservoir characterization and hydromechanical model calibration.

Seismic geomechanics (Task 1.1 and 1.2)

Feasibility studies

A feasibility/sensitivity study was carried out to calculate time-lapse seismic responses from a dynamic elastic reservoir model based on a North Sea deep reservoir undergoing large pressure changes. The time-lapse seismic travel-time shifts and time strains calculated from the modelled and processed synthetic datasets (i.e. pre-stack and post-stack data) were in a reasonable agreement with the observations, and indicated the feasibility of using a 1D strain rock physics transform and time-lapse seismic processing methodology. Estimated vertical travel-time shifts for the overburden and most of the reservoir were within ±1 ms of the observed values, indicating that the time-lapse technique is sufficiently accurate for predicting overburden velocity changes and hence geomechanical effects. Characterization of deeper structure below the overburden became less accurate, where more advanced time-lapse seismic processing and migration would be needed to handle the complex geometry and strong induced lateral velocity changes. Nevertheless, both migrated full-offset pre-stack and near-offset post-stack data imaged the general features of both the overburden and reservoir units. More importantly, results from the study indicated that integrated seismic and hydromechanical modelling can help constrain time-lapse uncertainty and hence reduce risk due to fluid extraction and injection (He et al. 2016).

Unlike reservoir models, which are commonly robustly calibrated with 4D seismic and fluids production data, geomechanical models are typically only loosely tied to available time-lapse seismic data, so we also undertook a study to assess what affect the uncertainty in typical input parameters has on modelled effective stress changes. Understanding how different modelling parameters affect stress change is extremely important for robust 4D seismic calibration. We used a Global Sensitivity Analysis technique on a sample of 1540 geomechanical models with parameters covering the range of typical of North Sea reservoirs, to assess the sensitivity of model inputs on the resultant vertical effective stress change. We successfully mapped each parameter to a sensitivity space and argue that by doing this we can reduce our initially considered ‘important’ inputs by 50%. Thus, reducing the number of calibration parameters and gaining a greater understanding of the model parameter sensitivity to stress and displacement allowed for a better understanding of what and how geomechanical models can be calibrated using time-lapse seismic data (Price et al. 2016a).

Anisotropy and fractures

Hydrocarbon production generally results in observable time-lapse changes within a compacting reservoir, but those physical property changes that lead to induced seismic effects, such as velocity increase, can be difficult to isolate uniquely. Thus, integrated hydro-mechanical simulation, stress-sensitive rock physics models and time-lapse seismic modelling workflows can be employed to study the effects of reservoir compaction on seismic velocity and seismic anisotropy. We studied the influence of reservoir compaction and compartmentalization on time-lapse reflection amplitude variation with offset (AVO) and with azimuth (AVOA). Specifically, the time-lapse AVO and AVOA responses were predicted for two reservoir models: a laterally homogeneous four-layer dipping model and a laterally heterogeneous graben. Seismic reflection coefficients for different offsets and azimuths were calculated for compressional (P–P) and converted shear (P–S) waves using an anisotropic ray-tracer as well as using approximate equations for AVO and AVOA. The simulations helped assess the feasibility of using time-lapse AVO and AVOA to evaluate induced stress anisotropy due to changes in the effective stress field. The results indicate that time-lapse AVO and AVOA analysis is a potential method for qualitatively and semi-quantitatively linking azimuthal anisotropy changes to pressure/stress change (He et al. 2015).

Fractures have a significant influence on the physical response of the subsurface. The presence of coherent fracture sets often leads to observable seismic anisotropy enabling seismic techniques to remotely locate and characterise fracture systems. We confirmed the general scale-dependence of seismic anisotropy and provided new results specific to shear-wave splitting (SWS). We found that SWS develops under conditions when the ratio of wavelength to fracture size (λS/d) is greater than 3, where Rayleigh scattering from coherent fractures leads to an effective anisotropy such that effective medium model (EMM) theory is qualitatively valid. When 1 <λS/d < 3 there is a transition from Rayleigh to Mie scattering, where no effective anisotropy develops and the SWS measurements are unstable. When λS/d < 1 we observed geometric scattering and began to see behaviour similar to transverse isotropy. We found that seismic anisotropy is more sensitive to fracture density than fracture compliance ratio. More importantly, we observed that the transition from scattering to an effective anisotropic regime occurs over a propagation distance between 1 and 2 wavelengths depending on the fracture density and compliance ratio. The existence of a transition zone means that inversion of seismic anisotropy parameters based on EMM will be fundamentally biased. More importantly, we observed that linear slip EMM commonly used in inverting fracture properties is inconsistent with our results and leads to errors of approximately 400% in fracture spacing (equivalent to fracture density) and 60% in fracture compliance. Although EMM representations can yield reliable estimates of fracture orientation and spatial location, our results showed that EMM representations will systematically fail in providing quantitatively accurate estimates of other physical fracture properties, such as fracture density and compliance. Thus more robust and accurate quantitative estimates of in situ fracture properties will require improvements to effective medium models as well as the incorporation of full-waveform inversion techniques (Yousef & Angus in review).

The presence of coherent fracture sets often leads to observable seismic scattering enabling seismic techniques to remotely locate and characterise fracture systems. We examined the widening effect of wavelets due to scattering within a fractured medium by using several different approaches. We used different methods including the RMS envelope analysis, shear-wave polarisation distortion, differential attenuation analysis and peak frequency shifting to assess the scattering behaviour of parametrised models in which the propagation direction is either normal or parallel to the fracture surfaces. The quantitative measures showed strong observable deviations for fractures of size in the order of, or greater than, the dominant seismic wavelength within the Mie and geometric scattering regimes for both propagation normal and parallel to fracture strike. The results suggest that strong scattering is symptomatic of fractures having size on the same order of the probing seismic wave (Yousef & Angus in review). Based on the anisotropy and scattering analysis, we assessed the feasibility of inverting SWS measurements to quantitatively estimate fracture strike and fracture density assuming an effective medium fracture model. The results of the full waveform synthetics indicate that the source frequency of the induced microseismicity (or seismic source) is crucial in extracting reliable fracture parameters due to the relationship between scale length of the probing seismic wave and the fracture heterogeneity (i.e., size). Although the SWS results themselves were diagnostic of fracturing, the fracture inversion allowed placing constraints on the physical properties of the fracture system. For real microseismic datasets, the range in magnitude of microseismicity (i.e., frequency content), spatial distribution and variable source mechanisms suggested that derivation of fracture properties from SWS measurements is feasible. For the single seismic source case and optimum receiver array geometry, the inversion for strike had average errors of between 11% and 25%, whereas for density average errors were between 65% and 80% for the single fracture set and 30% and 90% for the double fracture sets. Improvements on resolving strike can be made by including more microseismic sources in the inversion process. Furthermore, the improvements in resolving fracture density (or stiffness) can be achieved using a more advanced inversion approach such as anisotropic tomography in which the medium can divided into different domains (Yousef & Angus 2016).

Noise reduction

Noise is a persistent feature in seismic data and poses challenges in optimising seismic images and physical interpretation of the subsurface. We analysed passive seismic data from the Aquistore storage pilot project permanent seismic array to characterise, classify and model seismic noise. We performed noise analysis for a three-month subset of passive seismic data from the array and provided conclusive evidence that the noise field is not white, stationary, or Gaussian; characteristics commonly assumed in most conventional noise models. We introduced a novel noise modelling method that provides a significantly more accurate characterisation of real seismic noise compared with conventional methods, which is quantified using the Mann–Whitney–White statistical test. This method is based on a statistical covariance modelling approach created through the modelling of individual noise signals. The identification of individual noise signals, broadly classified as stationary, pseudo-stationary and non-stationary, provided a basis on which to build an appropriate spatial and temporal noise field model. Furthermore, we have developed a workflow to incorporate realistic noise models within synthetic seismic data sets providing an opportunity to test and analyse detection and imaging algorithms under realistic noise conditions (Birnie et al. 2016).

Noise is particularly troublesome for passive seismic monitoring where it commonly masks microseismic events. We proposed a statistics-driven noise suppression technique that whitens the noise through the calculation and removal of the noise covariance. Noise whitening was shown to reduce noise energy by a factor of 3.5 resulting in microseismic events being observed and imaged at lower signal-to-noise ratios than originally possible, whilst having a negligible effect on the seismic wavelet. The procedure was shown to be highly resistant to most changes in the noise properties and has the flexibility of being used as a stand-alone technique or as a first step before standard random noise attenuation methods (Birnie et al. 2016).

Rock Physics

Here we studied the uncertainty of various rock physics models in accounting for the non-linear relationship between effective stress and seismic velocity and assessed the associated error in velocity prediction. We used a collection of over 200 core measurements of ultrasonic velocity versus stress to constrain four end-member models: empirical, first principle, microstructural and a third-order elasticity model. We found all models provide a relatively good fit to the observed data, but some failed to accurately fit simultaneously both P- and S-wave data due to their not accounting for an observed stress dependent V_P/V_S ratio. Based on Bayesian statistical analysis we found that all model parameters are very well constrained. Using global probability distributions, we estimated parameter errors to be less than 5% for all models which propagate to errors in typical velocity change estimates to be of the order of ~10%. However, we noted third-order elasticity parameter errors of 10 to 15% with associated velocity errors on the order of ~100%. These large errors resulted from an under-determined inverse problem caused by the simplification of the model for an isotropic rock. We attempted to correlate model parameters with core porosity, but found porosity alone is not enough to constrain unknown model parameters within a suitable accuracy. The results of this study showed that errors in estimated velocity caused by model constraints are potentially well below errors introduced by natural recorded noise or noise introduced by 4D acquisition and signal processing.

Application to real case-studies (Task 1.3)

Snøhvit (Norway offshore)

A coupled fluid flow-geomechanical model for Snøhvit was built from data provided by Statoil, who also facilitated the model building process by providing direct feedback, clarification or updated data (i.e. depth converted horizons). The models were the focus of studies on the seismic response to pressure and fluid saturation changes.

Aquistore (Canada onshore)

The first post-CO₂-injection 3D time-lapse seismic survey at Aquistore was carried out following the injection of 36 kt of CO₂ within the selected reservoir units between 3170 m and 3370 m depth .

Limits on the detection of injected CO₂ in real noise conditions.

The objective was to assess the ability of time-lapse surface seismics to detect CO₂ in reservoir at ~3000 m depth after limited injection of CO₂. To simultaneously maximize repeatability and optimise subsurface imaging, particularly in the reservoir, a ‘4D-friendly simultaneous’ processing flow was used. This is identical to that used for processing the pre-injection data and is divided into pre-stack and post-stack elements. The repeatability between the baseline and post- CO₂-injection survey was excellent with global normalised root-mean squar8e (NRMS) values of 1.13 for the raw pre-stack data. The global NRMS value, tracked over the processing sequence, decreased with each processing step to a global NRMS value of ~0.10 after the final cross-equalization step was applied (Roach et al. 2016a).

Based on the processing and time-lapse analysis of the baseline and post-injection datasets and comparison with two previous pre-injection surveys (Roach et al. 2017), the following can be concluded:

• Repeatability between the seismic surveys was excellent because of the use of permanent buried geophones and buried dynamite sources.
• There is no apparent progressive degradation of the repeatability of the dynamite shots between successive surveys, though variations in data amplitudes have been observed.
• Of the three distinct units of the reservoir, the Deadwood reservoir unit registered a significant amplitude difference above the background and was previously identified, through fluid replacement modelling, as being the region most sensitive to changes in CO₂ content.
• The Deadwood anomaly corresponds to ~18 kt of CO₂, half of the injected volume, and the extent of the plume compares reasonably with the area calculated based on porosity and 100% saturation.
• In the other two injection zones, the Black Island sand and the lower Deadwood unit, no significant amplitude anomalies were observed.
• For the Black Island sand, repeatability is excellent so the seismic data points to the presence of a small amount of CO₂ in the zone. For the lower Deadwood unit the repeatability is not as good as the other two units and so the quantity of CO₂ is not detectable above the noise.
• The ability to compare the post-injection time-lapse results to pre-injection results allowed for better interpretation of the data.
• These results validate the value of the permanent array of buried geophones as well as the repeated survey parameters. Also, the availability of the pre-injection analysis is invaluable.

The background time-lapse signal-to-noise level at the Aquistore storage site was characterised through the use of a pre-injection time-lapse analysis performed on two sparse 3D seismic datasets acquired at the site. The analysis revealed a lowest global NRMS of 0.07, computed on time-windows across the whole dataset. NRMS levels above this would be taken to indicate the presence of CO₂ in the reservoir.

Gassmann fluid substitution and 3D seismic forward modelling was used to investigate the conditions under which injected CO₂ can be detected above the defined minimum noise level. Noise-free synthetic seismograms were generated for the baseline model and for monitor models with the replacement of brine by CO₂ at various saturations within the reservoir. Wave Unix (Xu 2012) was used to generate the synthetic surface reflection seismic data from an explosive surface p-wave source using. Two sets of models were used: (i) a thickness-CO₂ saturation model; and (ii) a single-zone model. The thickness-CO₂ saturation model included CO₂ distribution through varying the thickness of the CO₂ filled layer within the reservoir. The thickness of the layer was varied in 1m increments up to 100% of the zone thickness. For the single-zone model the reservoir is treated as single 193 m thick layer where CO₂ replacement is within the entire layer. Synthetic datasets were created for a range of CO₂ saturations of the reservoir layer.

To assess the detection limits for the injected CO₂ in the reservoir under the background noise level at the site, the noise traces from the baseline and pre-injection monitor surveys were added to each of the synthetic datasets. NRMS values from the CO₂ model synthetics with noise added were compared with NRMS values derived by comparison of the baseline and pre-injection monitor surveys within the same region. The CO₂ was defined as detectable when the NRMS from the models exceeded the NRMS value (0.07) from the pre-injection surveys. It should be noted that the NRMS of the synthetic dataset with added noise increases with increasing CO₂ saturation up to 20%, but it remains unchanged for further increases in saturation. Thus, with this metric, datasets with CO₂ saturations above 20% cannot be distinguished from each other.

We found that the time-lapse repeatability is excellent and provides the ability to monitor the CO₂ induced changes in the reservoir at the Aquistore storage site. Specifically, the analysis of the global NRMS for two different types of synthetic datasets with real noise added suggests that:

• For the thin layers, CO₂ is detectable in all zones under noise conditions provided that at least 5 m of the aquifer in each zone is saturated with 5% CO₂
• For 193 metre thick fully saturated layer it is possible to detect the CO₂ using a time-lapse analysis providing that the CO₂ saturation is above 20%.

This work is presented in Roach et al. 2016b. Note that the analysis considered only changes in reflectivity as a detectability criterion. Induced time-shifts are also an important time-lapse effect, particularly where CO₂ accumulates as a thicker layer, and can be very powerful detection tool.

Fracture characterisation using singular value decomposition (SVD)

The optimal basis formalism of Varela et al. (2007) was employed for the inversion of fracture density from synthetic AVOA data from a caprock within a geological setting with large impedance and anisotropy contrasts. The method is case-specific and is very dependent on knowledge of the in situ rock properties. In this case-study, the method is particularly advantageous because it is independent of the impedance contrasts of the layers. The model-based approach thus allows for direct quantification of the magnitude of anisotropy from an inversion of the p-wave amplitudes along the interface of interest within the high contrast geological setting.

The first stage saw the analysis of synthetic data using the optimal basis inversion to determine the fracture density of a set of synthetic data having different fracture densities characterising the magnitudes of anisotropy. The rock properties of the media are from well logs obtained from the Aquistore injection well. The effective horizontal transverse isotropy (HTI) media were defined by the slip-interface model of Liu et al. (2000) which explicitly includes the fracture density of the medium. The exact reflection coefficients at the interface between the upper isotropic medium and the HTI medium (fractured caprock) constituted the modelled data and contained a complete AVOA response of the reservoir. Atrak, an anisotropic ray tracer (Guest & Kendall 1993), was used to generate the synthetic surface reflection seismograms from which the amplitudes to be analysed were formulated (the amplitudes picked along the horizon between the caprock and layer immediately above it).

Two sets of synthetic data were evaluated: (i) a simple two-layer geometry which was modelled with the upper layer being isotropic described by the rock properties of the layer directly above the caprock and the lower layer being HTI with background rock properties of the caprock; and (ii) a 14-layer simplified geological model where the caprock was the sole HTI medium.

The second stage was to apply the optimal basis method to real data acquired at the Aquistore site. Three geological units were analysed: (i) the Winnipeg Icebox (the caprock) so as to determine any potential risks of leakage at the site and to establish the pre-injection state of the caprock; (ii) the Prairie Evaporite, located above the caprock, forming the regional seal and expected to be free of fractures and (iii) the Birdbear which is the shallowest horizon and is known to be fractured. The Prairie Evaporite and the Birdbear results form end-members and so provide a basis for assessing results from the caprock.

The final stage of this analysis was the determination of the direction of anisotropy for a complete characterisation of the fractures at the caprock. The Grechka & Tsvankin (1998) ellipse method was used.

The model-based optimal basis inversion method, combined with the use of an azimuthally dependent Shuey-type equation for extending the range of offsets for the inclusion the anisotropy signature, is suitable for determining the magnitude of anisotropy in the strong contrast deep reservoir environment at the Aquistore storage site. The crucial elements and/or possible limitations of this approach for the inversion of real data are:

• The ability to select the correct variations of in situ rock parameters in establishing the modelled reflectivity matrix in order to optimise the basis functions.
• The ability to extract true amplitudes in the real data as a function of offset/incident angles since the method weighs heavily on the assumption that modelled data is representative of the real environment.

A first look at the results from the real Aquistore data indicates different average magnitudes of anisotropy of the end members as expected — the Prairie Evaporite has an average crack density of 0.01 while the Birdbear has an average crack density of 0.02.

In Salah (Algeria onshore)

AVOA analysis

Elastic models from BP’s geomechanical simulations at the In Salah injection site (by Rob Bissell and provided by James Verdon) were observed to have negative elastic tensor components were observed where they should only be positive. James Verdon was then asked to provide the geomechanical outputs so that we could generate the elastic models from our codes. Unfortunately, using the geomechanical output (stress, Young’s modulus, pressure, etc) with our rock physics model transforms also generated negative values where they should only be positive. Attempts were made to track down the source of this problem but without success. However, we do not believe the source of the problem is our rock physics model as this issue has not been encountered with other geomechanical models.

Generating seismic synthetic data for the saturation models was not done because seismic acquisition parameters were not obtainable from BP.