We have carried out two seismic physical experiments to acquire wide-azimuth P-wave 3-D seismic data with a scaled down model (1:10 000) and scaled-up frequencies (10 000:1). Our aims are to verify the physical basis of using P-wave attributes for fracture detection, to understand the usage of these attributes and their merits, and to investigate the effects of acquisition geometry and structural variations on these attributes. The base model consists of a fractured layer sandwiched between two isotropic layers (Epoxylite). Inside the fractured layer there is a dome and a fault block for investigating the effects of structural variations. The two experiments were carried out using different acquisition geometries. The first experiment was conducted to maximize the data quality, with an offset-depth ratio of only 0.68 to the bottom of the fracture layer. For comparison, the second experiment was carried out to maximize the anisotropy effects, with the offset-depth ratio to the bottom of the fracture layer raised to 1.34. For each experiment, about 20 km2 of wide-azimuth 3-D data were acquired with a P-wave source. The physical modelling confirms that the P-wave attributes (traveltime, amplitude and velocity) exhibit azimuthal variations diagnostic of fracture-induced anisotropy. For the first experiment with noise-free data, the amplitude from the top of the fracture layer yields the best results that agree with the physical model parameters and free of the acquisition footprint. The results from other attributes (traveltime, velocity, AVO gradient) are either contaminated by the structural imprint, or by the acquisition footprint due to the lack of offset coverage. For the second experiment, despite the interferences from multiples and other coherent noise, the traveltime attributes yield the best results; both the acquisition footprint and the structural imprint are reduced due to the increased offset coverage. However, the results from the amplitudes are affected by the noise and are less reliable. Analysis of the two experiments reveals that the offset-depth ratio to the target is a key parameter for the success of the P-wave techniques. Smaller offset-depth coverage may only be applicable to amplitude attributes with high quality data; whilst large offset coverage makes it possible to use traveltime attributes. A reliable estimation from traveltime attributes requires an offset-depth ratio of 1.0 or more.
Fractures play an important part in reservoir development and enhanced oil recovery in tight formations of otherwise low permeability. The use of seismic anisotropy to detect natural fractures has been gradually gaining the acceptance of the hydrocarbon industry. The underlying physics for this technology comes from the equivalent medium theory for seismic wave propagation in fractured media, which has been intensively studied by a range of authors (e.g. Hudson 1981; Schoenberg & Douma 1988; Liu et al.2000a; Chapman 2003, amongst others). According to these theories, a medium containing vertically aligned fractures with scale length much less than the wavelength can be modelled by an equivalent azimuthally anisotropic medium for seismic wave propagation. Numerical modelling based on the equivalent medium theories reveals shear wave splitting and azimuthal variations in P-wave amplitude and traveltime as diagnostic features of the fractured medium.
Crampin (1981, 1985) was among the first in proposing the use of shear wave splitting for fracture characterization. Shear waves have been shown to be much more sensitive than P waves to detect fractures (Li 1997). There are two main reasons: first shear wave splitting can be measured from a single ray path, secondly the polarization of the fast shear wave and time-delay between the fast and slow waves are directly measurements of the fracture orientation and density. In contrast, P-wave measurements are made from a range of different ray paths in both the azimuth and offset domains (Ruger 1997; Li et al.2003). Furthermore, wide-azimuth 3-D P-wave data were not readily available in the 1980s. However, as the increasing use of wide-azimuth 3-D P-wave data in the late 1990s, coupled with the declining in shear wave acquisition, the use of P-wave seismic data for fracture detection became popular with the industry.
Over the past 10 yr, there has been a continual increase in the use of 3-D P-wave data for fracture characterization. Both numerical modelling and case histories of fracture detection using P-wave seismic have been the subject of intensive study (e.g. Lynn et al.1996; Liu et al.2000b; Smith & McGarrity 2001; Hall & Kendall 2003; Wang & Li 2006, amongst others). In comparison, the number of corresponding physical-modelling studies is much less (Luo & Evans 2004), whilst physical modelling studies of shear wave splitting were well documented in the early 1990s (e.g. Ebrom et al.1990; Brown et al.1991; Cheadle et al.1991; Slack et al.1991; Gregovic et al.1992). Here, we fill this gap by presenting physical modelling studies of fracture detection using large-scale 3-D P-wave seismic data. In addition to an examination of the underlying physics, this study also investigates the effects of acquisition geometries and compares the use of different P-wave seismic attributes and different analysis techniques on controlled experiments for fracture detection.
The Physical Models
As shown in Fig. 1, the base model consists of three horizontal layers. The first and third layers are constructed from the same material (epoxylite) and are believed to be isotropic. The second layer is constructed from a special industrial material and is believed to be azimuthally anisotropic. The material is composed of epoxy-bonded fibre sheets, simulating vertical fractures. The layer is highly anisotropic with about 20 per cent P- and S-wave anisotropy, and fracture density around 0.2 (Fig. 1). The measured elastic parameters and the corresponding anisotropic parameters are shown in Table 1, which exhibits some weak orthorhombic symmetry. There are also two built-in geological features inside the fracture layer. One is a dome, and the other is a fault block, consisting of a normal and a vertical fault (Fig. 1). The model is constructed with a scale of 1:10 000 for spatial dimensions and time measurements with a corresponding velocity scaling of 1:1.
Two models are derived from the base model for experiments simulating varying offset-depth ratios and acquisition geometries, as shown in Fig. 2. Model 1 (Fig. 2a) consists of a very deep-water layer of 1470 m on top of the base model in order to maximize the data quality, which ensures primary reflections from the base model free of multiple contaminations. The total thickness of the overburden above the fracture layer is 1962 m, and the maximum offset-depth ratio is about 0.9 to the top of the fracture layer, and about 0.7 to the bottom of the fracture layer. Model 2 consists of a thin water layer of 10 m and another isotropic layer of 430 m on top of the base model (Fig. 2b). Model 2 is designed to maximize the anisotropic effects. The total thickness of the overburden in Model 2 is reduced to 932 m, and the maximum offset-depth ratio is about 2.2 to the top of the fracture layer, and about 1.3 to the bottom of the fracture layer.
3-D data acquisition is conducted in a water tank. Fig. 3 shows the experimental set-up. The modelling system is shown Fig. 3(a), consisting of an ultrasonic pulse source and receiver system, analogue/digital converter and motor-driven positioning system. The maximum movement in the x, y and z directions are 230 × 230 × 100 cm, respectively, and the position error is less than 0.1 mm. The dominant pulse frequency is 230 kHz with a data bandwidth of 110–320 kHz. The source—receiver layout is designed to ensure a wide-azimuth coverage (Fig. 3b). Four shots are located in the centre of the spread, and are fired to 12 receiver lines with 768 channels. After firing four shots each time, the spread is moved to the next location. The receiver lines (Y-direction) are parallel to the fracture plane and the shot lines (X-direction) are perpendicular to the fracture plane (Figs 1 and 3). For Model 2, the same style of shooting was used but the acquisition parameters (e.g. source interval, receiver interval, etc.) are changed for increasing efficiency. Table 2 shows a detailed comparison of the acquisition parameters for Models 1 and 2.
For each experiment (or model), a total of about 20 km2 of 3-D data was acquired, and Fig. 4 shows a comparison of sample shot gathers from the two experiments. For Model 1, the data are of the highest quality, and all primary reflections are free from noise contamination. For Model 2, the primaries are contaminated with multiples and refracted arrivals, but the data are of reasonable quality. The top and bottom reflection events of the fractured layer are still clearly visible. Both data sets have a good coverage over azimuth (Fig. 5), and thus can be used for fracture analysis. The analysis is carried out in two stages: (1) initial data processing and (2) analysis of azimuthal anisotropy.
Initial Data Processing
For Model 1, the data are very clean and the processing is relatively straightforward. The flow includes: geometry loading, muting direct arrivals, CDP sorting, velocity analysis, NMO correction, stacking and post-stack migration. Good quality image volumes are obtained, and Fig. 6(a) shows a sample inline section. Both the dome and the fault block are clearly imaged. The image volume is used to guide traveltime and amplitude picking in pre-stack gathers.
For Model 2, a predictive model-based noise removal technique (Qian & Zhao 2003) was applied to the data before any other processing. Fig. 7 shows a comparison of a shot gather before and after noise removal. The data are reasonably clean after noise removal and the same processing flow as Model 1 is applied to the data afterwards. A sample inline section is shown in Fig. 6(b) for comparison. Due to the presence of noise, the section is not as good as Fig. 6(a) where the primaries are free from multiple interferences.
Azimuthal VariatiOns of P-Wave Attributes
The physical model consists of only a single set of fractures. Numerical modelling of wave propagation in such a medium predicts that the P-wave seismic attributes, such as traveltime, velocity and reflected wave amplitude will vary with azimuth, diagnostic of fracture-induced anisotropy (Liu et al.2000b; Hall & Kendall 2003). To verify this, we examine the azimuthal variations in P-wave attributes for both models. For Model 1, the maximum offset-depth ratio is 0.9 to the top of the fracture layer, and it is reduced to 0.7 to the bottom of the layer. In this case, we examine whether azimuthal variations in traveltime and velocity are observable. We then apply the same analysis for data acquired to Model 2 for comparison. The results are shown in Figs 8, 9 and 10, respectively.
Fig. 8(a) shows a super gather binned into 18 azimuthal gathers with 10° bin size. NMO correction is applied to these azimuthal gathers using a single velocity function aiming to flatten the event corresponding to the top interface of the fracture layer.
The top event of the fracture layer is properly flattened after NMO correction, whereas the bottom event shows azimuthal variation of residual moveout (Fig. 8a). For azimuths −50° and 40° and their adjacent gathers, the bottom event is also reasonably flattened (Fig. 8a). For azimuth 0° and its adjacent gathers, the bottom event shows under-correction (central panels of Fig. 8a). However, for azimuths −80° and 80°, as well as their adjacent gathers (left- and right-most panels, Fig. 8a), the bottom event is overcorrected. This indicates that the fracture normal is at azimuth 0° (X-axis), and the fracture strike is at azimuth 90° (Y-axis), since the wave propagates slower along the fracture normal than along the fracture strike. At an intermediate azimuth, the wave propagates with a velocity faster than at the fracture normal but slower than at the fracture strike. When an intermediate velocity for an intermediate azimuth is used to correct the moveout, it will overcorrect the events at the strike direction, but under-correct the event at the normal direction. This physical modelling result confirms these observations; it is more accurate and consistent compared with previous studies due to the improved data quality and experiment setup. However, since the maximum offset-depth ratio to the bottom of the fracture layer is only about 0.7 for this model (Table 2), there is a lack of offset coverage. As a result, the azimuthal variation in residual moveout is not fully developed.
Fig. 9(a) shows azimuthal offset-stacks for two offset ranges (0–800 m) and (800–2000 m) for three selected CDPs. For both offset ranges the azimuthal variation in traveltime is weak. Consequently the corresponding azimuthal velocity variation is also weak (Fig. 10a). There is no sufficient resolution in the velocity spectra to resolve the azimuthal variation, giving rise to a weak velocity ellipse as shown Fig. 10(a). The above confirms that azimuthal variations in traveltime and velocity can be observed for the data acquired for Model 1 which agrees with the equivalent medium theory. However, due to a small offset-depth ratio, the variation is not sufficiently developed. It is for this reason that a second experiment is conducted for Model 2 that has an increased offset-depth ratio.
For comparison, Figs 8(b), 9(b) and 10(b) show the azimuthal variations for Model 2 corresponding to Figs 8(a), 9(a) and 10(a) for Model 1. For Model 2, the offset-depth ratio to the top of the fractured layer is raised to 2.2, and that to the bottom of the layer is raised to 1.34 (Table 2). Compared with Fig. 8(a), the switch from under-correction to overcorrection is fully developed in Fig. 8(b) as a result of the increased offset-depth ratio, hence the increased offset coverage. As shown in Fig. 8(b), at azimuths −50° and 40°, the bottom event is properly aligned; at azimuth 0° and its adjacent gathers, the bottom event is under-corrected; but at azimuths −80° and 80°, as well as their adjacent gathers, the bottom event is overcorrected. The sinusoidal variation in the common-offset stack is also more profound and significant (Fig. 9b). Consequently, the azimuthal variation in the stacking velocities is also much clearer, giving rise to a strong velocity ellipse (Fig. 10b).
Numerical modelling also predicts that azimuthal variations of the P-wave attributes can be approximately described by an ellipse, as shown in Liu et al.(2000b) and Liu (2003). For velocity variation, the long axis of the ellipse indicates the fracture strike, and the ratio of the short to long axis is proportional to the fracture density. This can be clearly observed in Fig. 10(b) for Model 2. In comparison, the corresponding variation in Fig. 10(a) for Model 1 is much weaker due to a smaller offset-depth ratio.
To sum up, the above physical modelling confirms the numerical findings based on the equivalent medium theory. However, sufficient offset coverage is required to reveal these variations. The experiment shows that at least an offset-depth ratio of 1.0 is required to quantify these azimuthal variations reliably. Note that there exist a variety of equivalent medium theories for porous fractured medium (e.g. Hudson 1981; Sayers & Kachanov 1995; Liu et al.2000a, etc.). Although our numerical modelling is based on the equivalent medium theory of Hudson (1981) and Liu et al.(2000a), the findings here are not influenced by the choice of the corresponding theory. Liu et al.(2000a) show that the effects of the different types of theories on wave propagation are very small and negligible.
Extracting Fracture Parameters
As shown previously, the azimuthal variations of P-wave traveltime and velocity can be approximately described by an ellipse, and this is also true for interval traveltime as shown in Li (1999) and for amplitude and AVO gradient (Ruger 1997). For traveltime or interval traveltime, the short axis of the ellipse indicates the fracture strike, and for velocity, the fracture strike is indicated by the long axis. For the amplitude attributes, this depends on the impedance contrast. For a low-to-high impedance contrast, such as in the physical model studied here, the long axis of the ellipse indicates the fracture strike; for a high-to-low impedance contrast, the fracture strike is indicated by the short axis. In all cases, the ratio of the short to long axis is proportional to the fracture density. Therefore, we will use the fractional changes between the short and long axes to represent fracture density, that is, fracture density is equal to the difference between the long and short axes divided by the long axis [fracture density = (long-short)/long].
Two methods may be used to extract the fracture information from the P-wave attributes: full-azimuth surface fitting and narrow-azimuth stacking. The first method fits an elliptical surface to data from all available azimuths and offsets by a least-square fitting technique. The second method divides the data into a number of narrow-azimuth volumes; here we choose 18, with 10° azimuthal bins. Corresponding to these two methods, there are four principal seismic attributes: traveltime, amplitude, velocity and AVO gradient, which may be used to extract the fracture information. The surface fitting method is suitable for the traveltime and amplitude attributes. The narrow azimuth method is suitable for the velocity and AVO gradient attributes.
Both methods require the picking of traveltimes and amplitudes of the top and bottom of the target horizons through the pre-stack volume. Manual picking is difficult due to the workload and possible human picking errors. Thus, an automatic picker is usually employed. To ensure reliable traveltime and amplitude picking, the horizons are first manually picked from the post-stack volumes and then used as control points for pre-stack automatic picking. All traveltime and amplitude attributes are picked in this way.
The above methods and attributes analysis are applied to both data sets. We compare these results in order to evaluate the merits of different methods and different attributes, and assessing the effects of multiple interferences (noise), acquisition and structural variations on these attributes.
Fig. 11 shows the results of four attributes for comparison: top amplitude (amplitude of the top reflection event—interface 2 in Fig. 2) and its AVO gradient, interval traveltime and bottom traveltime (traveltime of the bottom reflection event-interface 3 in Fig. 2) of the fracture layer. Column (a) of Fig. 11 is for Model 1 and column (b) for Model 2. Except the gradient attribute, the other three attributes are analysed using the surface-fitting method; the AVO gradient is analysed using the narrow-azimuth stacking method. The result of the interval velocity is not shown since it resembles the traveltime. For each attribute in Fig. 11, we refer it to as Figs (11-a1), (11-b1), etc., indicating column (a) or (b), and row (1), etc.
For Model 1, the surface-fitting method applied to the top amplitude attribute gives the best results (Fig. 11-a1). The fracture strike is along the Y-axis, and the fracture density is about 0.2, with a very uniform distribution, indicating a single set of fractures. This agrees with the physical model. Also, the underlying structural features, as expected, do not affect the results. In contrast, the AVO gradient attribute by narrow-azimuth stacking does not give satisfactory results (Fig. 11-a2). The estimated fracture strike and fracture density are similar to those in Fig. (11-a1), but the distribution shows stripes parallel to the X-axis, indicating the effects of the acquisition footprint, enhanced by the binning and stacking process due to the lack of offset coverage.
The results from the interval traveltime by the surface-fitting method show a very clear structural imprint of the dome and the fault block (Fig. 11-a3). In areas outside these two zones, the estimated fracture parameters show very uniform distribution with the fracture strike along the Y-axis and an average fracture density of 0.2, agreeing with the physical parameters very well. The results from the bottom traveltime show similar features to the interval traveltime (Fig. 11-a4).
In contrast, for Model 2, the results from the amplitude attributes are affected by the presence of noise in the data, and are not reliable (Fig. 11-b1). The orientation estimation varies more significantly compared with those in Fig. (11-a1) for Model 1, and the estimated fracture density is scattered due to noise interference. The binning and stacking process further enhances the noise effects, as shown in Fig. (11-b2).
Instead, the surface-fitting method applied to the bottom taveltime attribute gives the best results (Fig. 11-b4), with the fracture strike along the Y-axis, and with a very uniform distribution of fracture density. Also the effects of the structural imprint are compensated for by the increased offset-depth ratio to the fractured target. This is because the magnitude of the anisotropy variation is greater than the structural variation once the offset-depth ratio reaches 2.0. The interval traveltime gives out similar results but slightly affected by the presence of the fault (Fig. 11-b3).
Analysis and Findings
We analyse and compare the results in order to draw some useful guidelines related to the use of P-wave seismic for fracture detection.
Choice of attributes
The amplitude is the most sensitive attribute. For clean data, even when the offset-depth ratio to a target is very small (x/z = 0.9), fracture parameters can still be determined accurately from the amplitude attributes (Fig. 11-a1, Model 1). However, the presence of a small amount of noise will significantly distort the results (Fig. 11-b1, Model 2) even with sufficient offset coverage (x/z = 2.0). Therefore, for the use of the amplitude attribute, it is important to reduce the noise and preserve the amplitude than increase the offset coverage.
The use of traveltime attributes required sufficient offset coverage to allow the azimuthal variation sufficiently developed. For insufficient coverage (offset-depth less than 1.0), the attributes will be heavily influenced by the structural imprint and also by the acquisition footprint (Figs 11-a3 and 11-a4, Model 1), and these effects are substantially reduced when the offset coverage are doubled (Figs 11-b3 and 11-b4, model 2).
Choice of methods
The surface fitting method to all offsets and all azimuths is preferred than narrow-azimuth stacking. The later enhances the acquisition footprint when the offset-depth ratio is not sufficient large (Fig. 11-a2, Model 1), and also enhances the noise response even the offsets are doubled (Fig. 11-b2, Model 2).
Effects of acquisition
The offset and azimuthal coverage is critical for the success of using P waves for fracture detection. The offset-depth ratio to the target is a good indicator. The experiment for Model 1 reveals that the offset-depth ratio should be at least larger than 1.0. Only in a noise-free environment, may there be sufficient sensitivities from the amplitude attribute to resolve the fracture parameter for offset-depth ratio (x/z) approaching 0.9 (Fig. 11-a1, Model 1). Large offset coverage makes it possible to use the traveltime attributes. A reliable estimation from traveltime attributes required offset-depth ratio of 1.0 or more (Figs 11-b3 and 11-b4).
Effects of structural variation
The structural variation will leave a significant imprint on the estimated results, particularly if the offset-depth ratio is not sufficient large (Figs 11-a3 and 11-a4, Model 1). However, these effects may be compensated for by increasing the offset coverage and by the use of traveltime attributes. When the magnitude of azimuthal variation due to anisotropy exceeds that due to structural variations, the structural imprint will be significant reduced (Figs 11-b3 and 11-b4, Model 2).
We have carried out two physical experiments and performed a detailed analysis of two 20 km2 3-D data sets acquired from the experiments. Azimuthal variations of P-wave attributes are observed, confirming the numerical modelling results based on the equivalent medium theory. Two methods (full-azimuth/full-offset and narrow-azimuth/full-offset) have been used and four seismic attributes analysed for fracture parameter estimation. The results from the narrow-azimuth stacking method appear to be influenced by the acquisition footprint and the noise, but those from full-azimuth and full-offset surface fitting agree with the physical parameters.
The offset-depth ratio to the target is a key parameter for the success of the P-wave techniques. It affects the choices of attributes and choice of processing methods. Smaller offset-depth coverage may only be applicable to amplitude attributes with high quality data; whilst large offset coverage makes it possible to use traveltime attributes, which is less sensitive to the presence of multiple interference, reducing the effects of acquisition footprint as well as the structural imprint. A reliable estimation from traveltime attributes requires offset-depth ratio of 1.0 or more.
We are grateful to CNPC for providing the funds to conduct the physical experiments, and for permission to show the results. We thank Wang Shoudong for processing the data and Mark Chapman and Enru Liu for comments on manuscript. The work is supported under a collaboration agreement between the CNPC Geophysical Key Laboratory at the China University of Petroleum and the Edinburgh Anisotropy Project (EAP) at the British Geological Survey (BGS), and is published with the approval of all project partners. This work is also supported by the National Scientific Foundation of China under contract of 40574056, and the Ministry of Education under the New Century Elite Scientist Support Scheme (NCET-04-0106). The work is presented with the permission of the Executive Director of the British Geological Survey (NERC).