Angular Dispersion: Seismic Anisotropy, Thomsen Parameters, and Fracture Detection
Angular dispersion is the variation of seismic wave velocity as a function of propagation direction, encompassing both azimuth (compass bearing within the horizontal plane) and dip (angle from vertical or from horizontal). In a perfectly homogeneous, isotropic medium, seismic waves travel at the same speed in every direction and angular dispersion is zero; any measurable deviation from directional uniformity constitutes anisotropy, and angular dispersion quantifies the magnitude and geometry of that deviation. The phenomenon arises from three principal geological mechanisms in sedimentary basins: the preferred orientation of clay mineral platelets and mica flakes in fine-grained rocks such as shale and mudstone, which creates faster horizontal than vertical P-wave propagation (vertical transverse isotropy, or VTI); the presence of aligned vertical or sub-vertical open fractures that transmit seismic energy preferentially parallel to the fracture strike (horizontal transverse isotropy, or HTI); and the coexistence of both mechanisms in the same rock, producing orthorhombic anisotropy with nine independent elastic constants. Angular dispersion is quantified using the Thomsen parameters: epsilon (ε) describes the fractional difference between horizontal and vertical P-wave velocities, gamma (γ) describes the equivalent contrast for shear waves polarized horizontally versus vertically, and delta (δ) governs the angular dependence of P-wave velocity between zero and 90 degrees from vertical, controlling the ellipticity of the wavefront at intermediate angles. Typical shale formations in the WCSB display epsilon values of 0.05 to 0.25, gamma values of 0.05 to 0.20, and delta values ranging from minus 0.05 to plus 0.15. Angular dispersion affects every aspect of seismic data quality and interpretation: it distorts depth conversion when using vertical velocities measured by sonic logs to predict lateral propagation speeds, biases normal-moveout velocity analysis toward intermediate values between the fast and slow axes, and creates apparent structural relief on migrated images if the migration velocity field fails to account for the directional variation in wavefront propagation speed.
Key Takeaways
- VTI anisotropy from clay alignment in shale: Vertically transverse isotropy is the most common form of seismic anisotropy in sedimentary basins and arises when clay mineral platelets and micas align preferentially in the horizontal plane during compaction, creating a medium that is elastically equivalent to a stack of infinitely thin isotropic layers. In VTI media, P-wave velocity is fastest in the horizontal direction (along layering) and slowest in the vertical direction (across layering), with the angular dependence of velocity controlled by the Thomsen parameters epsilon and delta. The Duvernay shale formation of central Alberta, a primary unconventional reservoir with P-wave vertical velocity of approximately 5,200 m/s and horizontal velocity of approximately 5,700 m/s, displays epsilon values of 0.10 to 0.18 as measured from walkaway VSP surveys and multicomponent core ultrasonic testing. Failing to account for this 10 to 18 percent velocity difference when converting seismic travel times to depth causes systematic depth errors of 50 to 180 metres at Duvernay targets near 3,500 to 4,000 metres depth, mislocating well landing points and compromising the value of expensive horizontal wells.
- HTI anisotropy from aligned vertical fractures: Horizontal transverse isotropy arises when a set of parallel, sub-vertical open fractures is embedded in an otherwise isotropic or VTI host rock. Seismic waves traveling parallel to the fracture planes propagate faster than waves traveling perpendicular to them, because waves crossing the fractures must repeatedly cross compliant, fluid-filled or gas-filled crack boundaries that slow propagation. The result is azimuthal anisotropy: P-wave velocity and amplitude vary with the azimuth of the source-receiver direction relative to the fracture strike. In the Devonian carbonate formations of the southern Alberta Foothills, natural fracture sets oriented northeast-southwest can produce P-wave azimuthal velocity anisotropy of 2 to 8 percent, sufficient to create clearly measurable amplitude and travel-time differences between inline and crossline traces on 3D seismic data. Mapping HTI anisotropy magnitude and fast-axis azimuth provides a direct estimate of fracture strike that can be integrated with image log data from horizontal wells to design perforation cluster spacing and hydraulic fracture stage spacing.
- Shear wave splitting as a fracture indicator: When a shear wave enters an HTI medium, it splits into two orthogonally polarized components: the fast shear wave (S1), polarized parallel to the fracture strike, and the slow shear wave (S2), polarized perpendicular to the fracture strike. The time delay between the arrivals of S1 and S2 at a receiver accumulates as the waves travel through the fractured interval and is proportional to the fracture intensity (crack density) and the thickness of the fractured layer. Shear wave splitting is the most sensitive seismic indicator of fracture-induced anisotropy and can detect crack densities as low as 0.02 (one crack per 50 crack-radii-cubed of rock volume). Multicomponent seismic data acquired with full-vector receivers and multicomponent source vibrators are required to record shear wave splitting; conventional P-wave surveys can detect only the azimuthal P-wave velocity anisotropy that is a secondary consequence of fractures. In Duvernay exploration, where natural fractures are a key uncertainty in well performance, multicomponent 3D surveys with shear wave splitting analysis cost approximately CAD 2 to 4 million more than a conventional P-wave program of the same size but provide fracture intensity maps that can improve the success rate of horizontal well placement by 15 to 25 percent.
- AVAZ analysis for fracture azimuth and intensity: Amplitude variation with azimuth (AVAZ) analysis uses the azimuthal dependence of P-wave reflection amplitudes on 3D seismic data to estimate the orientation and intensity of sub-seismic fracture sets without requiring multicomponent data. At a given angle of incidence, the P-wave reflection amplitude from a fractured layer varies sinusoidally with source-receiver azimuth: it is largest when the survey azimuth is perpendicular to the fracture strike (where the HTI anisotropy most strongly affects the AVO gradient) and smallest when parallel to the fracture strike. By fitting an ellipse to the azimuth-amplitude relationship at each CMP location in a 3D dataset, interpreters extract the fast axis azimuth and the anisotropy magnitude, both of which are proxies for fracture orientation and intensity. AVAZ analysis requires dense azimuthal sampling: the 3D acquisition must be designed with multiple source and receiver azimuths at each CMP, achieved through wide-azimuth or full-azimuth shooting geometries with source-receiver azimuth bins of 22.5 to 45 degrees providing adequate sampling for stable ellipse fitting.
- Depth imaging and Thomsen parameter estimation: Kirchhoff pre-stack depth migration and reverse-time migration both require an accurate anisotropic velocity model to correctly focus seismic energy and position reflectors at true depth. In VTI media, even a small positive delta of 0.05 causes the NMO velocity measured from surface seismic to exceed the true vertical velocity by approximately 5 percent at offsets equal to the depth, producing apparent uplift of 50 metres on a 1,000-metre-deep reflector if isotropy is assumed. Estimating Thomsen parameters for the migration velocity model requires integrating well-log P-wave and S-wave velocity anisotropy measurements from ultrasonic core plugs or cross-dipole sonic logs, walkaway VSP data that sample the velocity field at intermediate propagation angles, and surface seismic velocity analysis at multiple offsets and azimuths. Full-waveform inversion with anisotropic parameterization can update the Thomsen parameter fields simultaneously with the interval velocity model, but requires broadband seismic data from 2 to 80 Hz and is computationally intensive, with run times of days to weeks on high-performance computing clusters for WCSB-scale 3D datasets of 200 to 500 km2.
Measuring and Correcting for Angular Dispersion in WCSB Seismic Programs
Measuring angular dispersion from surface seismic begins with sorting 3D CMP gathers into azimuth sectors, typically 22.5-degree wide bins covering the full 180-degree azimuth range from 0 to 180 degrees (northwest-southeast to northeast-southwest in most WCSB programs). Within each azimuth sector, the normal-moveout velocity is picked on key reflection events as a function of offset, producing a velocity field that varies with azimuth if anisotropy is present. The azimuthal variation of NMO velocity follows an ellipse, with the fast azimuth corresponding to the direction of maximum horizontal velocity and the slow azimuth to the minimum. Fitting this ellipse to the azimuthal velocity picks yields the fast-azimuth direction (interpreted as the fracture or stress orientation) and the anisotropy magnitude (the ratio of fast to slow NMO velocity), which for typical HTI formations in the WCSB Devonian ranges from 1.01 to 1.08 (1 to 8 percent).
Correcting seismic data for VTI angular dispersion in depth migration requires building a model with four parameters per layer: vertical P-wave velocity (V0), epsilon, delta, and the vertical P-to-S velocity ratio. Vertical velocity is constrained by sonic logs in wells within the 3D footprint. Epsilon and delta are estimated from the azimuthally averaged NMO velocity field: if the short-offset NMO velocity approximates V0 / sqrt(1 + 2 x delta) and the long-offset moveout approximates V0 x sqrt(1 + 2 x epsilon), then differencing the two velocity estimates at a given horizon gives access to the epsilon-delta contrast. Walkaway VSP surveys, which fire a surface source at multiple offsets and azimuths while recording on a borehole geophone array, provide the most direct constraint on Thomsen parameters because they sample intermediate propagation angles that surface seismic cannot reach with practical offset ranges.
The practical consequence of ignoring angular dispersion is most severe in the Foothills, where dipping Paleozoic carbonate targets beneath thrust sheets require accurate anisotropic velocity models for both imaging and depth conversion. A depth error of 100 metres in a Foothills Devonian prospect translates directly into a well cost error: each 100 metres of additional horizontal displacement from the planned well trajectory requires approximately CAD 80,000 in additional drilling time plus the risk of missing a thin target entirely if the casing point is set at the wrong depth. Foothills Devonian wells typically cost CAD 8 to 14 million, making the CAD 400,000 to CAD 800,000 incremental cost of an anisotropic tomography velocity update a sound investment against the risk of a depth miss.
In wide-azimuth Duvernay 3D programs, AVAZ analysis is routinely applied to map fracture intensity variations across the reservoir section. A 150 km2 wide-azimuth 3D program acquired with 8 azimuth sectors and maximum offsets of 5,500 metres provides sufficient azimuthal sampling for AVAZ gradient extraction at the Duvernay level near 3,700 metres depth. The AVAZ-derived fracture intensity map is integrated with isotropic AVO gradient maps to produce a composite fairway map, highlighting areas of both gas saturation (from AVO) and natural fracture density (from AVAZ) that represent the highest-value horizontal well landing zones. In the Willesden Green Duvernay fairway of west-central Alberta, AVAZ programs have identified northeast-trending fracture corridors of 1 to 3 km width where natural fracture intensity doubles the estimated stimulated reservoir volume per horizontal well stage, increasing expected 18-month cumulative production from approximately 200,000 to 350,000 barrels of oil equivalent compared to unstimulated fracture positions.
Fast Facts
Thomsen epsilon values for Duvernay shale range from 0.10 to 0.22, meaning horizontal P-wave velocity exceeds vertical velocity by 10 to 22 percent and introduces depth errors of 50 to 180 metres at typical Duvernay depths of 3,500 to 4,500 metres if isotropic migration is used. Shear wave splitting time delays in moderately fractured Devonian carbonates of the WCSB typically range from 5 to 30 milliseconds per kilometre of travel path through the fractured interval, corresponding to crack densities of 0.03 to 0.12. P-wave azimuthal velocity anisotropy detectable by AVAZ analysis requires a minimum azimuthal NMO velocity contrast of approximately 0.5 percent, equivalent to a 15 m/s velocity difference at V0 = 3,000 m/s, which is achievable with signal-to-noise ratios greater than 5 on pre-stack gathers at the target horizon. Wide-azimuth 3D seismic programs in the Duvernay play cost 20 to 35 percent more than conventional orthogonal single-azimuth programs of the same area because of the additional shot lines, receiver templates, and extended acquisition time required to fill all azimuth-offset bins uniformly.