Constraints on fluid flow pathways from shear wave splitting in and around an active fluid-escape structure: Scanner Pockmark, North Sea

SUMMARY

Vertical fluid-escape structures observed in seismic reflection data represent an important class of potentially active fluid flow pathways. An understanding of the mechanism of fluid flow in these types of structures is needed to assess the risk of natural gas venting from potential subsurface carbon dioxide storage operations. The Scanner Pockmark Complex is a 22 m deep, 900 × 450 m seabed depression in the North Sea, which actively vents methane, and is underlain by a seismic chimney structure with horizontal dimensions of ∼300 × 600 m. Gas accumulation is evidenced by the presence of bright reflectors at the top of this seismic chimney, at a depth of ∼50 m below the seabed. Here, we analyse seismic anisotropy in these shallow sediments using shear wave splitting observed on ocean bottom seismographs (OBS). Anisotropy varies spatially, with a strength of ∼1–4 per cent, on several OBS located in and around the pockmark complex. By correlating these observations with calculated subsurface P- and S-wave velocities, we show that there is anisotropy present throughout the sediments through which the chimney passes, which are interpreted as relating to syn- and post-depositional glaciomarine processes. However, within the chimney itself the orientation of the fast direction is different to that outside the chimney and the degree of anisotropy is lower. We attribute this difference as indicating that the anisotropy observed within the chimney is associated with the formation and continued presence of the gas migration system, which overprints the background depositional anisotropy.

Fracture and flow, Gas and hydrate systems, Controlled source seismology, Seismic anisotropy

1 INTRODUCTION

1.1 Chimney structures

Chimney structures (Hustoft et al. 2010; Moss & Cartwright 2010a) represent an important class of seismically imaged fluid-escape features. They are observed in seismic reflection data as vertical to subvertical anomalies with circular or elliptical planforms, displaying seismic blanking and/or discontinuous or chaotic reflections (e.g. Løseth et al. 2011). Chimneys have been observed in sedimentary basins and on rifted margins globally (e.g. Cartwright et al. 2007; Gay et al. 2007; Moss & Cartwright 2010a,b; Løseth et al. 2011; Plaza-Faverola et al. 2017), including extensively in the North Sea (Hovland & Sommerville 1985; Cole et al. 2000; Bünz et al. 2003; Karstens & Berndt 2015). For example, analysis of 3-D seismic reflection volumes in the South Viking Graben, North Sea (an area of 2850 km²; Karstens & Berndt 2015) identified 46 large-scale (∼100–1000-m-wide) chimneys located within the shallowest 1000 m of the overburden.

Large-scale chimney structures (∼100–1000 m wide) have been suggested as regions where there is elevated permeability relative to the normal ‘background’ permeability of the host sediment (Cartwright et al. 2007; Cartwright & Santamarina 2015), for example by a series of interconnected subvertical or radial fractures, which allow the vertical flow of gas in the shallow subsurface (Bull et al. 2018). However, there is an important distinction between the appearance of ‘seismic’ chimneys and true ‘physical’ chimney structures. A seismic chimney is generated as a result of the large energy partitioning of a gas accumulation in the subsurface, which produces a zone of acoustic blanking below a gas layer giving the appearance of a vertically oriented chimney or pipe. However, this observation does not identify or demonstrate the actual gas migration pathway in this region (the ‘physical’ chimney), which may not match the geometry of the apparent seismic manifestation. For example, sealing and lateral spreading processes at the top of the gas accumulation may cause the blanking zone to appear to cover a greater area than any primary migration pathway, resulting in an apparent seismic chimney larger than the physical chimney, or even where there is no physical chimney present at all.

Time-lapse seismic experiments conducted prior to, during, and after subsurface gas release in the QICS (Quantifying and Monitoring Potential Ecosystem Impacts of Geological Carbon Storage, located in Ardmucknish Bay, Oban, Scotland; e.g. Cevatoglu et al. 2015; Taylor et al. 2015; Räss et al. 2018) and STEMM-CCS experiments (Strategies for Environmental Monitoring of Marine Carbon Capture and Storage, located in the North Sea; Flohr et al. 2021; Roche et al. 2021) have been used to observe the formation mechanism and structure of physical chimneys in the shallow overburden (Cevatoglu et al. 2015; Bull et al. 2018; Roche et al. 2021). In settings characterized by shallow depths, low effective stress, and fine-grained and poorly consolidated sediments, chimney formation is expected to occur primarily by the generation of new fractures or flow pathways, with the first phase of formation being the hydraulic fracturing of low permeability sediments, or re-activation of pre-existing fractures due to high fluid overpressure (Arntsen et al. 2007; Cartwright et al. 2007; Løseth et al. 2009; Fauria & Rempel 2011). With increasing effective stress conditions, found at greater depths due to the increased overburden thickness, upward percolation of fluids may also occur through capillary flow processes (Cathles et al. 2010).

Fractures in sedimentary settings play an important role in the properties of subsurface reservoirs, as they enhance porosity and permeability, or conversely may contribute to reservoir compartmentalization. Thus, fracture orientations, densities and sizes are of interest to the understanding of subsurface fluid flow, both for natural fluid escape systems and potential subsurface CO₂ storage operations (e.g. Robinson et al. 2021), as well as hydrocarbon reservoirs (e.g. Odling et al. 1999; Aydin 2000; Bratton et al. 2006) and geothermal settings (e.g. Jafari & Babadagli 2011; Ghassemi 2012). Attribute analysis of seismic image volumes, which includes techniques such as coherence analysis (e.g. Bahorich & Farmer 1995), are able to detect larger fractures. However, these techniques are unable to image smaller fractures below the spatial resolution of the seismic image. Analysis of seismic anisotropy provides an alternative approach to determine fracture properties associated with chimney structures. This technique uses the directional dependence of transmitted seismic wave speeds. Several theories have been developed to relate the elastic response of fractured rocks to fracture distribution and properties (e.g. Hudson 1981; Thomsen 1995; Liu et al. 2000; Chapman 2009; Jakobsen & Chapman 2009; Jin et al. 2018).

1.2 Anisotropy from shear wave splitting

The measurement of seismic anisotropy using shear wave splitting (SWS) is an established technique for determining the presence and orientation of fracture networks (e.g. Hudson 1981; Crampin 1985; Schoenberg & Douma 1988; Simmons 2009). Thomsen (1986) provides a framework in which several anisotropic parameters are defined using observations from seismic data and assembled to form a stiffness tensor for the anisotropic symmetry system in question. Using this approach the characteristic patterns produced in seismic data by different types of anisotropy can be determined. The two dominant end-member anisotropic systems are: horizontal transverse isotropy (HTI), which is associated with vertically aligned features, such as vertical cracks, and; vertical transverse isotropy (VTI), which dominates where the primary linear structure is horizontally aligned, for example in stratigraphic layering.

In the study of active source seismic anisotropy in relatively shallow settings, SWS analysis typically uses P-to-S converted seismic waves (e.g. Li et al. 1988; Haacke & Westbrook 2006; Exley et al. 2010; Tsuji et al. 2011). When an incident polarized shear wave enters an anisotropic medium, either as an S-wave from an isotropic medium or during conversion from the P-wave at the boundary, it is split into two separate shear components (Fig. 1a). One component is polarized parallel to the primary linear fabric of the medium (e.g. cracks/fractures, crystal preferred orientation, stress orientation), while the other is, in the simplest cases, oriented perpendicular to this. The fabric-parallel component has a faster wave speed than the fabric-perpendicular component (Fig. 1a, labelled 1 and 2). By analysing the patterns recorded by seismographs, the orientations of the polarizations of the two shear waves, and, hence, the symmetry axis, can be determined. The delay time between the arrival of the fast and slow waves is used to determine the intensity of anisotropy, which is related to the fracture size and/or density (e.g. Crampin 1985; Mueller 1992; Li 1997). Where multiple successive layers are anisotropic, the eventual arrivals at the instrument will represent a summation of the effects of each layer. If the orientation of the symmetry plane is different between successive anisotropic layers then repolarization will occur, further splitting the incident fast and slow waves into additional fast and slow pairs for each (Fig. 1a, label 3), and resulting in more complex composite patterns of SWS.

Figure 1.

Shear wave splitting (SWS) in an HTI medium. (a) Schematic diagram of SWS in an HTI medium. P-wave paths from the source (star, located ∼10 m below sea surface) to the OBS located on the seafloor (orange triangle) are shown as solid black lines. S-wave legs are shown as dashed lines. Polarized shear waves entering or generated by conversion at the base of anisotropic layers are split into two orthogonal components, oriented in the fast and slow (red and blue arrows) anisotropic directions. Where multiple layers are anisotropic (e.g. labels 1 and 2), the eventual arrivals at the instrument will represent a summation/superposition of the effects (label 3). (b) Plan view definition of components and geometry used for analysis. Black arrows show orientations of X and Y components. Red and blue arrows are directions of S1 and S2 components respectively. In this study we set S1 as the component which has a rotation angle in the quadrant 0–90° clockwise from north. Purple and green arrows show radial and transverse directions for an example shot, shown as a cross. Shots are located around a constant offset circle around the instrument (orange triangle). (c)–(k) Synthetic seismic record sections showing SWS patterns for a simple 2-layer case of vertically oriented HTI, following co-ordinate definitions in (b). Symmetry plane azimuths are 70/160° for the upper anisotropic layer and 30/120° for the lower anisotropic layer. Time = 0 s is the direct arrival. Only the direct arrival and P–S converted arrivals are shown. Azimuths are measured clockwise from north. (c) OBS X and (d) Y components. Horizontal geophones are oriented N–S and E–W respectively. Rotation of X and Y components into (e) fast (S1) and (f) slow (S2) component orientations. (g) Slow component with time shift applied to remove delay between fast and slow components (event 4). (h) Radial component generated from initial X and Y components shown in (c) and (d). Event labelled 6 shows the time-varying nature of the radial component arrival from an anisotropic layer. (i) Radial component following layer stripping correction. Arrival labelled 6 is now flattened. (j) Transverse component conjugate to (h). Events labelled 7 show the characteristic set of four polarity reversals and amplitude nulls at 90° intervals. (k) Transverse component conjugate to (i), following layer stripping. Energy has been removed from the arrivals (event 7). Following layer stripping to correct for the anisotropy in the upper layer, blurred edges of the 4 × 90° polarity pairs in (j) become clear and aligned at 30/120° in (k), blue dashed lines are plotted for reference.

In this study, we are concerned with understanding the role of fracturing in vertical fluid flow associated with subsurface chimney structures. Therefore, we look for evidence of HTI-type anisotropy. The characteristic patterns of HTI SWS are best observed by transformation of the two perpendicular horizontal geophone components (Fig. 1c and d) into an alternative co-ordinate system comprising: (1) a radial component (Fig. 1h) where the direction of energy transmission from each shot is in the direction from the shot to the receiver and (2) a transverse component (Fig. 1j) where the energy transmission is perpendicular to this. The idealized pattern produced by HTI symmetry systems comprises: (1) on the radial component, an azimuthally continuous arrival which shows a cos2θtraveltime variation with shot-receiver azimuth with a pair of peaks and troughs each separated by 180° (Fig. 1h, event 6) and (2) on the transverse component, two pairs of 90°-wide opposing polarity arrivals with their centres located at 180° from each other, separated by amplitude nulls at 90° intervals (Fig. 1j event 7). The orientation of the polarity reversals and amplitude nulls gives the orientations of the two perpendicular symmetry planes in the HTI anisotropic system. Here, and throughout this article, we use the term ‘event’ to mean a group of arrivals on a record section.

If the observed anisotropy results from fracturing, rock physics may be used to provide insights into the fracture properties. A range of theories have been developed to describe the elastic response of fractured rocks (e.g. Hudson 1981; Thomsen 1995). While these generally agree for dry rock, they differ considerably where fluids and fluid flow between cracks and pores are present (Liu et al. 2000; Amalokwu et al. 2014, 2015a,b, Jin et al. 2018), as may be the case for fluid flow through chimneys. In addition, uncertainties arise because a medium containing a small number of large fractures will generate the same response as a medium containing a larger number of smaller cracks (e.g. Maultzsch et al. 2003).

SWS analysis of active source seismic data recorded on ocean bottom seismometers (OBS) has previously been used to study marine slope stability at a number of locations, including on the West Svalbard continental slope (Haacke & Westbrook 2006), to identify vertical fluid migration pathways within the gas hydrate bearing sediments in the Storegga slide offshore Norway (Exley et al. 2010), to understand the state of stress in the seismogenic subduction zone of the Nankai Trough, Japan (Kimura et al. 2021), and in the context of conducting accurate depth imaging of subsurface reservoirs (Sil et al. 2010).

1.3 Aims

In this study, we test the hypothesis that the structure of the primary fluid pathways within a subsurface chimney structure can be detected and measured using seismic anisotropy. We apply this technique to the study of a large-scale active fluid-escape system in the Central North Sea, the Scanner Pockmark Complex. We compare observations of seismic anisotropy made within, around, and away from the pockmark and chimney, to determine the effect of the fluid flow conduit on the ‘background’ geology. This is done by identifying the orientations of the fast and slow directions and the magnitude of SWS. We correlate observations of SWS with the subsurface stratigraphic structure by combining observations from seismic reflection imaging and modelling of shallow subsurface P- and S-wave velocity structure. The key aim of this study is to further understand the potential structural controls on fluid-escape at the Scanner Pockmark Complex, and, by extension, other actively venting chimney structures.

2 STUDY AREA

2.1 Scanner Pockmark

In this study, we investigate the chimney structure observed below the Scanner Pockmark Complex (SPC), located in the Witch Ground Basin of the North Sea (Fig. 2). Pockmarks are formed during erosive fluidization, or blow-out events, when fluid overpressure in the subsurface reaches a critical level, leading to rapid upwards migration (Hovland et al. 2002). The SPC is a compound structure comprising two large pockmarks; the West Scanner and East Scanner pockmarks, each of which are >75 m wide, >250 m long and >15 m deep. Strong acoustic backscatter within the water column provides evidence for present-day active venting at the Scanner pockmarks, enabling a calculated present-day gas flux in the range 1.6–2.7 × 10⁶ kg yr^–1 (272–456 L min^–1; Li et al. 2020). The presence of methane-derived authigenic carbonates which arise due to oxidation of the escaping gas (Judd et al. 1994; Judd & Hovland 2009) supports the persistence of this venting over time.

Figure 2.

Scanner Pockmark Complex study region and experiment geometry. a) Location of North Sea licence box 15/25, containing the SPC. Dashed black line demarcates boundary between United Kingdom (to W) and Norway (to E) exclusive economic zones (200 nm). (b) Ship-acquired swath bathymetry of the study area, around Scanner Pockmark Complex and additional class 1 pockmarks (Scotia, Challenger) to the north. Primary pockmark complexes are labelled in bold, with East and West Scanner Pockmarks labelled in italics. Class 2 pockmarks are visible as small features across the whole region, with a NNW–SSE orientation. Dashed black box indicates region shown in (c). (c) Scanner Pockmark OBS experiment geometry. Inverted red triangles show OBSs used for SWS analysis and discussed in this paper. Inverted white triangles show other OBS locations, with red outline indicating instrument was used in shallow velocity analysis (see Appendix  B, note OBS 1 was not used in velocity analysis and has black outline). Grey filled triangles indicate OBSs which fail initial QC and tilt analysis tests and are excluded from further analysis. OBS numbers are labelled. Dark grey dashed lines show ship track during GI source acquisition, used for investigation of SWS. Light grey dotted–dashed lines show ship track during Duraspark surface sparker acquisition, used for subsurface imaging and velocity model construction (see Section 5). Solid black lines indicate locations of Duraspark and SBP profiles shown in Fig. 3.

The region surrounding the SPC contains a large number (>1500 over an area of 225 km²) of smaller pockmarks (Figs 2b and c), which are interpreted as forming due to localized pressure changes and sediment dewatering (Böttner et al. 2019), as well as several further large pockmarks, the Scotia, Challenger and Alkor Pockmarks (Gafeira & Long 2015).

2.2 Geological history, seismostratigraphic framework and regional stress regime

The Scanner Pockmark Complex is hosted within the uppermost part of the ∼600 m-thick Quaternary sediments of the Witch Ground Basin (Fig. 3), which are described in detail by Stoker et al. (2011), Böttner et al. (2019) and Callow et al. (2021). The principal stratigraphic subdivisions, listed from bottom to top, with their corresponding Marine Isotope Stages (MIS), are:

1. The Aberdeen Ground Formation (AbG; MIS 100–13), which comprises layered sands, silts and clay-rich sediments and has a well-layered seismic character.

2. The Ling Bank Formation (LB; MIS 12–10), which is interpreted to represent a sub-facies characteristic of a glacial tunnel valley (e.g. Kluiving et al. 2003; Graham et al. 2007; Ottesen et al. 2014), with valley widths up to 1.5 km and depths of 100 m (Callow et al. 2021), and consisting of: a lower subunit comprising coarse sands and gravels, and displaying a chaotic seismic character; an intermediate clay-rich subunit which is seismically homogenous, and; an upper subunit comprising coarse sands and characterized by higher amplitude reflections.

3. The Coal Pit Formation (CP, MIS 6–3), which comprises muddy sands interpreted to represent glacial tills.

4. The Last Glacial Maximum Formation (LGM; MIS 3–2), which comprises silty-sandy clays with rare pebbles.

5. The Witch Ground Formation (WG; MIS 2–1), which is composed of silty clay sediments and is itself subdivided into two units:

• (5.1) The Fladen Member, which we refer to here as the Lower Witch Ground (LWG), and which displays an interbedded seismic character.

• (5.2) The Witch Member, or Upper Witch Ground (UWG), which displays a uniform seismic character.

Figure 3.

Subsurface reflection images of geological structure beneath Scanner Pockmark and reference site. Profiles are taken along black lines in Fig. 2(c). (a–c) Hull mounted sub-bottom profiler. The sub-bottom profiler data shows reflection strength of arrival, black = reflections, white = zero amplitude. (d–f) Duraspark surface sparker seismic data. Seismic data has + and - polarities (black/white), zero amplitude is in grey. Inverted white triangles show locations of OBSs on the seafloor along these profiles, labelled with instrument number above. Vertical white dashed lines show outline of seismic chimney in the Aberdeen Ground Fm.. Dashed white ellipses at top of LB (0.28 s TWTT) show locations of high amplitude reflections interpreted as gas accumulation, which is predominantly in the LB and breaches upwards into the CP (Callow et al. 2021). Dashed white ellipse at 0.25 s TWTT in (e) shows laterally restricted high amplitude reflection which may represent migratory gas in the physical chimney just below the pockmark. Abbreviations for stratigraphic units: UWG, Upper Witch Ground Fm.; LWG, Lower Witch Ground Fm.; WG, Witch Ground Fm.; LGM, Last Glacial Maximum Fm.; CP, Coal Pit Fm.; LB, Ling Bank Fm. and AbG, Aberdeen Ground Fm..

The Scanner Pockmark depression erodes down to the base of the Witch Ground Fm.. Beneath the Quaternary stratigraphy are the Palaeogene and Neogene Hordaland and Nordland Groups, which are composed of low-permeability claystone (Judd et al. 1994).

In the Witch Ground Graben, the present-day regional maximum horizontal stress, σ1, is orientated NW–SE (Zanella & Coward 2003), and the horizontal stress exceeds the vertical stress. The minimum horizontal stress direction, σ3, is ∼54° (Evans & Brereton 1990), and therefore, if present, tensional fractures would be expected to form perpendicular to this, that is at ∼140–150°. Mapping of subsurface faulting and glaciomarine structures and orientations was performed by Callow et al. (2021).

2.3 Location of chimney and gas

High-amplitude zones in the seismic reflection surrounding the SPC are interpreted as gas-saturated sediment layers, and are observed at three discrete stratigraphic levels (Callow et al. 2021): (1) the Crenulate Reflector (CR) located at the base of the AbG (i.e. the base of the Quaternary); (2) within subunit 2 (following the definition of Callow et al. 2021) of the LB, where it breaches through vertical fluid pathways into the base of the CP and (3) at the base of the WG, immediately below the pockmark.

Analysis of seismic reflection data (Figs 3d and e; Callow et al. 2021) shows a vertical seismic chimney beneath the SPC. This feature appears to have an elliptical planform, displaying a diameter of ∼600 m in the N–S direction and ∼300 m in the E–W direction, which is consistent with the orientation of the long axis of the Scanner Pockmark planforms. The seismic chimney extends through the AbG. The accumulation of gas at the top of the chimney in the LB can be seen as a region of high amplitude seismic reflections above the blanking zone. Below the seismic chimney, the CR is laterally continuous across the region (Callow et al. 2021), indicating that any effects of gas-blanking from the accumulation at the LB are reduced by this point.

Mapping of the RMS amplitude below the base of WG in 3-D seismic data allows the lateral distribution of gas in the LB to be mapped (Callow et al. 2021). Gas is present immediately below both of the Scanner pockmarks, and spreads ∼400–500 m from the centre of West Scanner in a quadrant between NW and NE. A further accumulation extends toward the NW, initially with a narrow width (<200 m across) but widening with increasing distance from the pockmark. Gas is also inferred to be present at depth beneath all of the large pockmarks in the area surrounding the SPC: the Scotia, Challenger and Alkor pockmarks (Callow et al. 2021), although it is not limited to solely being beneath the pockmarks. The locations of the pockmarks represent an interaction between gas presence at depth and subsurface geological drivers.

3 DATA ACQUISITION

RRS James Cook cruise JC152 conducted a seismic experiment over the Scanner and Challenger Pockmark Complexes. The seismic data were recorded by an array of 25 ocean bottom seismographs (OBS; Fig. 2c), each equipped with a hydrophone (hereafter notated as H) and three-component (X, Y, Z) geophone package, and recording at a sampling rate of 4 kHz, corresponding to a sample interval of 0.25 ms. Eighteen instruments were deployed around the Scanner Pockmark Complex, with OBSs 1 and 2 located within the West Scanner Pockmark and OBS 3 within the East Scanner Pockmark. A further seven OBSs were deployed at a reference site located ∼1 km southeast of the Scanner Pockmark Complex, in an area where subsurface reflection imaging showed there to be no evidence for subsurface gas accumulation or migration.

The experiment used four different seismic sources, GI and Bolt airguns and Duraspark and Squid 2000 surface sparkers, to provide a broad bandwidth and permit the investigation of frequency-dependence of anisotropy (Fig. 2; Bull 2017). To achieve maximal azimuthal coverage, which is necessary for determining the directionality of anisotropy, profiles were acquired at multiple orientations through the OBS array (Fig. 2c). In this study we analyse the recordings of the GI airgun source, which was acquired with a 420 ci (2 × 105/105 ci) array working in harmonic mode at a pressure of 2000 psi. In addition, we use recordings of signals generated by a Duraspark surface sparker operated at 2000 J. Arrivals from each of the seismic sources were recorded during a single phase of OBS deployment. Signals produced by the GI airguns and surface sparkers were also recorded on two towed multichannel streamers, and used in the development of the integrated stratigraphy beneath the SPC (Callow et al. 2021). The RRS James Cook was also equipped with a hull-mounted Kongsberg SBP120 sub-bottom profiler, which recorded data with a 2.8–6 kHz frequency bandwidth and 4.4 kHz central frequency, and a Kongsberg EM710 multibeam sonar system.

4 METHODS

Throughout the description of the data processing below, we use the following nomenclature to define the different OBS geophone components (Fig. 1b). X, Y and Z refer to the original instrument geophone component orientations, which are together generalized hereafter as G, where X and Y are the horizontal components and Z is vertical. G^T is used to indicate tilt-corrected components. G’ components are G^T components which have undergone a horizontal rotation about the true vertical Z axis, where S1 and S2 represent a specific case of X’ and Y’ where these are oriented in the fast and slow anisotropic directions, and N and E represent a specific case where the components are oriented approximately north and east. The radial (R) and transverse (T) components are oriented in and perpendicular to the source–receiver directions, with the angles used for rotation defined such that R is oriented in the same relative direction for all instruments.

4.1 Data processing 1: quality control to rotation of geophone components

The first stage of data processing involves assignment of the survey geometry, data quality control (QC) and transformation of the horizontal geophone components into R and T components.

4.1.1 Shot and receiver locations

The first step in the processing sequence was to determine the true location of each OBS on the sea bed, as the initial co-ordinates only correspond to the drop locations at the sea surface. The actual OBS co-ordinates were estimated using a simple grid search algorithm to determine the locations that minimized the difference between the observed and modelled direct wave arrival times. For each OBS a solution was found that yielded an RMS traveltime difference between the modelled and observed arrival times of less than 1 ms for all source types. This corresponds to a location error of <1.5 m.

4.1.2 Data QC

Initial QC of the OBS records comprised tests for horizontal vector fidelity and frequency fidelity. Horizontal vector fidelity is tested by plotting the azimuth vs absolute amplitude of the X and Y component direct arrivals, which should show a |cos 2θ| pattern with peaks of approximately equal height. Frequency fidelity is determined by calculating frequency spectra for shot lines passing close to the instrument. The horizontal vector fidelity test is failed by OBS 3, while OBSs 4 and 17 show issues with their frequency spectra for the GI gun source. These instruments are excluded from further analysis. Channel polarity assignment is necessary to resolve ambiguity in the assignment of the polarities of the two horizontal geophone components. These were checked by hodogram analysis of first motions associated with straight shot lines, and are assigned to be consistent throughout the data set. During this stage of QC we also noted a high-amplitude oscillatory signal with a frequency of ∼9 Hz present on all OBSs.

During QC we noted that the X geophone component on several instruments appeared to show a delay to the peak of the direct arrival when compared to the Y and Z components (Fig. A1). This was also accompanied by a longer wavelet, indicating a lower frequency content, confirmed by frequency spectrum analysis (Fig. A1f). There is a positive correlation between the magnitude of the delay and the difference in frequency content. These differences appear to be a result of an instrument response. The effect of the observed behaviours on the data appears to be that of a filter, which we hereafter denote as the X geophone functional filter. To solve for the X geophone functional filter effect, we treat it as two separate components, a static shift, which can be measured from the data, and a zero-phase filter, the applications of which result in an improvement to the R:T power ratio (see Appendix  A for further information).

4.1.3 Geophone package tilt analysis

We performed a geophone package tilt analysis following the approach used by Exley et al. (2010). For each sample (0.25 ms) in a 3 ms window (encompassing the first half-cycle of the direct arrival) around each shot within an offset range of up to 500 m from the OBS, the shot-receiver azimuth, the direct arrival polarization azimuth and the difference between these two quantities are calculated. The standard deviation of the differences between the direct arrival polarization and shot-receiver azimuth are then calculated, and iterative 3-D rotations of the geophone components are performed to minimize this misfit using the rotation matrix

$$\begin{equation*} \left( {\begin{array}{@{}*{1}{c}@{}} {{N^T}}\\ {{E^T}}\\ {{Z^T}} \end{array}} \right)= \left( {\begin{array}{@{}*{3}{c}@{}} {\sin \alpha }&{ - \frac{{\cos \alpha \cos \beta }}{{\sin \alpha }}}&{ - \frac{{\cos \alpha \cos b}}{{\sin \alpha }}}\\ 0&{\frac{{\cos \gamma }}{{\sin \alpha }}}&{ - \frac{{\cos \beta }}{{\sin \alpha }}}\\ {\cos \alpha }&{\cos \beta }&{\cos \gamma } \end{array}} \right)\left( {\begin{array}{@{}*{1}{c}@{}} N\\ E\\ Z \end{array}} \right) \end{equation*}$$

α and β are the angles from the vertical of the N and E components (with 90° corresponding to horizontal) and γ is the angle of Z from the vertical calculated using α and β. The rotation is performed iteratively over a range of test values for α and β of ±10° from the horizontal, at intervals of 0.25°. Full details of geophone tilt analysis are provided in Appendix  A.

4.1.4 Rotation to radial-transverse

Transformation of the X^T and Y^T geophone components to R and T was used to identify the characteristic patterns of SWS associated with HTI anisotropy (Figs 1h and j). This process was achieved by finding the rotation angle for each shot which minimizes the energy of the direct arrival on the T component, which should theoretically be zero. This process thus maximizes the power ratio R:T (Figs A3f,g and A4f,g). The minimization calculation was performed on a trace-by-trace basis, resulting in a rotation angle being calculated for each shot (Figs A3d,e and A4d,e). As with the tilt analysis described above, the successful implementation of this process is highly dependent on the proper alignment of the input X^T, Y^T direct arrival wavelets. Discussion of the effects of the various corrections may be found in Appendix  A.

4.1.5 Summary—QC and rotation

We find that for the Scanner Pockmark seismic experiment the typical range of geophone package tilts, given as the angle of Z from the vertical, γ, is in the range ∼1.5–7.5°. This is consistent with the low seabed relief in the survey area. For these small values of γ we observe relatively little difference in the rotation angles or the R:T power ratio found using the R–T rotation calculation.

In theory, for two orthogonally oriented horizontal geophone components a single rotation angle should govern the transformation to R and T components. However, this is not always the case in practice, primarily due to residual errors in survey geometry. To produce the final R and T components we followed the approach of Exley et al. (2010) and Haacke & Westbrook (2006) and fitted a sinusoid of the form |$a + b\cos ( {c\theta - d} )$| to the trace-by-trace angles calculated by the power minimization routine. In this equation a represents the mean orientation azimuth of the X geophone component, b is the magnitude of the distribution around this mean angle, c is the period of the variability, which is expected to be ∼2 and d is a phase shift component. This approach has several advantages in that (1) it allows more reliable rotation of traces at offsets where the geometry is less well defined, (2) it permits flexibility in the application of various corrections (for example, if it is desired to not apply the filter component in removing the effects of the X geophone functional filter) and (3) it fits the trend in the data rather than the noise (Exley et al. 2010).

4.2 Data processing 2: flattening, binning, stacking and visualization of arrivals

4.2.1 Moveout and stacking

To correctly characterize azimuthal variations, the effects of variations in source–receiver distance must first be removed from the data—a process that we call ‘flattening’. In our first study of anisotropy with this data set (Bayrakci et al. 2021), only the very shallowest parts of the OBS data records were targeted. Flattening of the arrivals was performed using a simple geometry-based time correction based on the seismic velocity of water (1490 m s^–1) and the relocated shot-receiver geometry. This approach can be used when the upgoing S-wave leg is effectively vertical, since the converted P–S arrivals will be approximately parallel to the direct arrival. This condition holds for the shallowest part of the subsurface in settings where Vp/Vs is very large. However, for deeper conversion points, as Vp/Vs reduces (principally due to a faster increase in Vs relative to Vp) the S-wave leg is further from the vertical, and so a different flattening approach is required. Normal moveout (NMO) corrections can be used to apply variable flattening with depth.

However, since Vp > Vs, the ray path taken by the P–S converted arrivals is highly asymmetric (Fig. 1a), with the S-wave conversion occurring at a location close to directly beneath the instrument. As a result, the standard hyperbolic NMO equation does not fully hold, resulting in overcorrection at greater offsets. This non-hyperbolicity can be accounted for by the addition of a quartic term to the standard quadratic hyperbolic NMO definition. However, the coefficients of these terms must then be determined. Alternatively, the effects of non-hyperbolic moveout can be mitigated by restricting the velocity analysis to limited offset ranges. Here we take this latter approach. P–S moveout velocity functions were obtained for overlapping ∼200 m-wide offset ranges (50–200, 100–300, 200–400 m, etc.) for the radial components for each instrument, with the relevant offset range trends then applied to flatten arrivals (Fig. 4). We observed that the NMO velocity decreases from water values (∼1490 m s^–1) to values ≤1000 m s^–1, generally stabilizing at ∼0.4 s after the direct arrival. The non-hyperbolicity results in larger values of NMO velocity for the longer offset windows.

Figure 4.

Radial component semblance panels for a) OBS 1, and b) OBS 12. Semblance is calculated in offset windows of 50–200 m, 100–300 m and 200–400 m, which approximately correspond to the different offset ranges of square shot profiles around instruments (Fig. 2c). Data are filtered with a 20–30–60–100 Hz zero-phase bandpass filter. Dashed black lines show the picked P–S semblance trend in each case. Key characteristics of the semblance spectra are annotated. Note, times are shifted such that the seabed is located at 0 s. High semblance values indicate velocity-time pairings which produce good flattening of the arrivals.

4.2.2 Shot selection, binning and stacking

We took sections of shot profiles forming approximate squares to generate ‘walkaround’ shot profiles. The offset range of the included shots lies within as small a range and are distributed as symmetrically as possible to minimize the effects of differential normal moveout. Due to the acquisition geometry used in the experiment the offset ranges included in stacking are typically 100–200 m. At different offsets different time intervals of the records are covered by arrivals which are not related to P–S conversion processes, for example, by the direct arrival multiple. This is an issue in this study due to the shallow water depth and target of imaging. While shot locations exist at a range of offsets in this study, offset-depth balancing must also be considered since shorter offset arrivals would be expected to show better imaging potential for the shallowest layers. As the shot location and the conversion depth vary, the location of the conversion point moves and so the sampling location may also vary.

Shots were azimuthally sorted into fixed-width bins and then stacked to improve signal-to-noise ratio. Smaller bin widths resulted in a smoother azimuthal coverage and are more suitable for longer offset arrivals. For shorter-offset shots a larger binning interval is required to avoid empty bins. A bin width of 9° was selected as it represents a good balance between improving signal-to-noise ratio and achieving a full, smooth azimuthal coverage around the instrument. The azimuthal resolution is limited by the binning interval used, and depends on both the data quality and the acquisition geometry (e.g. shot and line spacing).

4.2.3 Layer stripping

Raw estimates of SWS give only depth-averaged estimates of anisotropy. Therefore, to determine the depth variation of both fracture orientation and density a layer stripping approach is required (Haacke et al. 2009; Exley et al. 2010), which recursively compensates and removes the anisotropy measured in each layer. To demonstrate this process, synthetic seismograms were generated using ANRAY (Gajewski & Pšenčík 1987), for a model comprising two anisotropic layers. Traveltimes were calculated through the anisotropic model subsurface by ray tracing, then seismograms were generated using a 200 Hz wavelet, which was subsequently bandpass filtered to match the frequencies at which signals are observed in the data set used in this study. The model has a water layer of 150 m thickness with Vp = 1480 m s^–1, an upper layer of thickness 15 m with Vp = 1640 m s^–1, and a lower layer of thickness 20 m with Vp = 1740 m s^–1. Vs values for each layer are calculated using the relationship of Castagna et al. (1985) for mudrocks. Below the lower anisotropic layer is an isotropic half-space. Both anisotropic layers have an anisotropy of 2.5 per cent. In Figs 1(c)–(k), the arrivals at 0 s are the direct arrival at the OBS, while the arrivals at ∼0.18 and ∼0.27 s correspond to the converted P–S arrivals from the first (label 1 in Fig. 1a) and second (labels 2 and 3 in Fig. 1a) layers, respectively.

Layer stripping is performed by first rotating the instrument X and Y horizontal geophone components (Figs 1c and d) in order to align the components with the fast and slow (S1 and S2 respectively) orientations for the shallowest layer (Figs 1e and f). These orientations are found by inspection of the Rand T components. This is followed by the application of a static time shift to align the target arrivals in the S1 and S2 components (Fig. 1g, event 4). Since the sign of the delay on the corresponding R component arrival is not always clear, we choose to use the symmetry plane orientation determined from the T component in the range 0–90° to perform this rotation from X, Y to S1, S2. We then determine which of S1 and S2 is the fast and slow component by analysing the sign of the required time shift. Following application of the time shift, the rotation was reversed from S1, S2 to X, Y. These components were then re-rotated to R and T to analyse the effects of removing the anisotropy on successively deeper layers (Figs 1i and k). When removing the shallowest observed anisotropy on a given instrument, the layer stripping correction flattens the target arrival on the R component (Fig. 1i, event 6), with the sinusoidal-type delay being transferred up the record to the previously flat direct arrival. On the T component, the amplitude of the anisotropic target arrivals is reduced, and energy is also transferred upwards to the direct arrival (Fig. 1j, event 7).

For the example synthetic model shown in Fig. 1, the orientation of the symmetry planes of the lower of the two anisotropic layers is rotated relative to the upper, being aligned 30/120° rather than 70/160°. The effect of this difference in alignment can be seen in the lower pairs of arrivals on the S1 and S2 components (Figs 1e–g, event 5), where the edges of the 180°-wide polarity pairs are blurred compared to the upper set (event 4). This can also be seen on Fig. 1(j), event 8, which shows blurring at the edges of the four 90°-wide sets of polarity reversals, compared to event 7 where the sharp edges are clearly visible and define the symmetry planes at 70/160°. The effect of applying the layer stripping correction for the upper anisotropic layer can be seen in Fig. 1(k), event 8, where the lower set of arrivals now do show sharp energy nulls at 30/120°, the orientation of the symmetry planes for this lower layer.

The time delay between the fast and slow arriving shear components may either be accumulated uniformly between the time of an observation of anisotropy and the next shallowest corrected event, or over a limited vertical range within this interval. This impacts the percentage anisotropy which would be calculated for a layer. If layer boundaries that cannot be resolved in the data are missed then the net result is that estimates for percentage anisotropy, particularly for the deeper subsurface layers, represent lower bound values.

The layer stripping approach described here holds only for HTI anisotropy, since the rotation from X, Y to S1, S2 is a 2-D rotation. It is possible to generate synthetic seismograms with more complex organizations of symmetry planes (Gajewski & Pšenčík 1987), for example by the inclusion of an inclination from the vertical (known as tilted transverse anisotropy, or TTI). These could then be used to identify characteristic patterns of more complex systems in the OBS R and T components. However, in order to layer strip tilted anisotropic cases an Alford (1986) rotation approach, which needs an additional orthogonally polarised source, is required (e.g. Winterstein & Meadows 1991; Tsvankin et al. 2010; Li et al. 2019).

5 CONSTRUCTION OF SUBSURFACE VELOCITY MODELS FOR STRATIGRAPHIC CORRELATION

To relate observations of converted waves (and splitting) in the OBS data to P-wave reflection data it is necessary to generate sufficiently finely depth-resolved models of subsurface P- and S-wave velocity. Velocity analysis was performed using a semblance approach applied to the hydrophone component recordings of the Duraspark surface sparker source. A sparker source was chosen for this analysis since it offers a higher resolution for the shallowest layers (which are the region of interest) and gives clear early-arriving near-offset (up to 500 m) waveforms. Instruments were selected for the analysis to give a range of settings and provide simple solutions. Therefore, instruments and shot lines associated with subsurface dipping structures at depth (Fig. 3) and instruments located above shallow gas accumulations were excluded. A total of nine instruments were used in velocity modelling (Fig. 2).

Moveout (flattening) velocities were picked using a combination of single shot-profiles and whole instrument gathers, with a single composite moveout trend determined for each instrument (Figs B1a and B2a). These were then converted to a layered P-wave interval velocity model using the Dix (1955) equation and converted from velocity–time to velocity–depth (Figs B1b and B2b). Forward modelling was performed with horizontally stratified blocks of constant velocity (no vertical velocity gradients). We chose to utilize a stepped model, that is one with only sharp velocity contrasts between layers and no vertical velocity gradients, since sharper/larger velocity contrasts are required to generate stronger reflections. Therefore, this type of model should represent a good fit to the primary reflectivity that is observed in the streamer data (Fig. 3).

Each P-wave velocity model was then tested by forward ray-tracing to generate predicted traveltimes for each of the reflection layers in the model. These predicted traveltimes were compared to picked arrival times from the hydrophone records and the models were iteratively improved if required (Figs B1c,d and B2c,d). To test the sensitivity of the velocity modelling the models were perturbed and fits recalculated. Further details of this process are provided in Appendix  B.

The resulting P-wave velocity–depth structures (Fig. 5) were converted to S-wave velocity using the empirical relationship of Castagna et al. (1985) for mudrocks:

$$\begin{equation*} {\rm{Vp}}\left( {{\rm{km }}\,\, {{\rm{ s}}^{-1}}} \right) = 1.16{\rm{Vs }}\left( {{\rm{km}} \,\,{{\rm{ s}}^{-1}}} \right) + 1.36 \end{equation*}$$

Figure 5.

Velocity structure at the Scanner Pockmark and reference site. (a) P-wave velocity and (b) S-wave velocity structures, as determined through analysis in Appendix  B. Black lines are profiles from OBSs located around pockmark, grey lines are OBSs located at reference site. Solid red line is the mean value calculated from all 9 profiles shown. Dashed red line is the range of values around the mean calculated from model sensitivity analysis, as discussed in Appendix  B. Light grey shading indicates total range of velocity values plus the calculated errors from sensitivity analysis. Stratigraphic units are labelled following the abbreviations in Section 2.2.

The calculated S-wave velocities at the seabed (Fig. 5) are in the range ∼60–110 m s^–1. These values are consistent with laboratory S-wave velocity measurements from samples collected at Scanner pockmark, which show a range of ∼55–100 m s^–1 for the uppermost few meters (Bayrakci et al. 2021), and with shear wave velocities derived from other locations in the North Sea (e.g. Armstrong et al. 2020).

To predict P–S converted wave arrival times and compare against the arrivals observed on the OBS records, the individual P- and S-wave velocity–depth profiles were averaged, allowing Vp and Vs to be determined for each stratigraphic unit. Two-way traveltime picks were then made beneath each OBS, to account for lateral variations in depth to different horizons arising as a result of, for example, regional dips and erosional unconformities, and then converted to depth using the P-wave velocity. The traveltime for the P–S arrival at zero offset was then calculated using the P- and S-wave velocities for the vertical downward and upward ray path legs respectively, from each picked interface beneath, assuming only the simplest conversion mode of a downgoing P-wave and an upgoing S-wave. Due to the averaging procedure applied to velocities, the most significant source of error in predicting P–S times arises from uncertainty in the S-wave velocity, which itself depends on the P-wave velocity uncertainty (discussed in Appendix  B) and the chosen relationship between Vp and Vs.

6 OBSERVATIONS

We present evidence for SWS on several OBSs located within the pockmark, adjacent to the pockmark and at the nearby reference site where a seismic chimney structure and gas venting are not present. Instruments included in the analysis are those which display good evidence for P- to S-wave converted arrivals and for HTI-type SWS patterns, such as those shown for an idealized, synthetic case in Figs 1(c)–(k).

For each instrument, we analyse the system from the shallowest layers down, applying layer stripping calculations to remove the effects of shallow anisotropy where required. For each instrument discussed we present a series of figures that show the evidence for the P–S converted arrivals and SWS present in the radial, transverse, and separated fast and slow oriented components. The vertical axis in each figure is in P–S time, which is a composite two-way traveltime for the converted wave energy. Azimuths are measured clockwise from north. For each figure, features of interest are labelled with numbered circles which correspond to the arrival events discussed in the text below. For the first instrument discussed, OBS 19, we also compare the observations to the idealized synthetic case in Figs 1(c)–(k).

When calculating the anisotropy, strictly it is necessary to consider two aspects of the propagation pathways of the P- and S-wave phases. First, the percentage anisotropy as a rock physical property depends on the orientation of the energy propagation relative to the causative structure of the anisotropy. Due to the high Vp:Vs ratio (Fig. 5) and relatively short offsets involved in this study, we assume that the upgoing S-wave legs are approximately vertical, and, therefore, parallel to the along-fracture direction if the anisotropy is pure HTI-type. Secondly, it is necessary to separate the traveltime accumulated as a downgoing P-wave and an upgoing converted S-wave, with the traveltime delay due to anisotropy occurring on the latter upwards leg. Because Vp >> Vs in this setting (Fig. 5) we can assume that the S-wave traveltime is approximately equal to the difference between the P–S and direct arrival traveltime, particularly for the shallower layers where the path length is the shortest. The effect of this assumption is that the percentage of anisotropy may be slightly underestimated. Since the ray paths are oriented approximately parallel to the along-fracture direction for HTI-type anisotropy (e.g. Rathore et al. 1995), in the following discussion we define percentage anisotropy as the ratio of the S1 to S2 traveltime difference to the S-wave traveltime through the layer of interest.

6.1 Reference site

To characterize the ‘background’ anisotropy across the study area, we first analyse evidence for anisotropy at the reference site located ∼1.5 km SE of the pockmark complex (Fig. 2c). By studying any patterns of anisotropy observed at this location we seek to understand regional anisotropy, over which any anisotropy associated with the gas migration system at the pockmark may be superimposed. Subsurface reflection images show no seismic chimney structures or significant gas accumulations in this location (Fig. 3f). Two instruments, OBS 19 and OBS 20, show good evidence of SWS in the Quaternary sedimentary sequence at this location.

6.1.1 OBS 19

On OBS 19, following the direct arrival multiple, which is seen as the upward curving arrivals at 0.16–0.18 s, we see characteristic HTI-type 90°-spaced amplitude nulls on the T component at ∼0.18–0.22 s, with reversals in the polarities of the adjacent energy highs either side (Fig. 6b, event 2, comparable to Fig. 1j, event 7). The orientation of the amplitude nulls/polarity reversals, and, hence, the anisotropic symmetry planes, is ∼30/120°. Over the same time interval we also observe a corresponding set of arrivals on the R component that are azimuthally continuous around the instrument, though they show some variation in their arrival time (Fig. 6a, event 1, comparable to Fig. 1h event 6). Flattened R component arrivals in this time interval could not be related to the direct arrival multiple due to the significantly lower NMO velocity used (e.g. Fig. 4). Above 0.16 s we do not see evidence for characteristic HTI SWS patterns. Stratigraphic correlation using the P- and S-wave velocity structures derived in Section 5 places the origin of these arrivals close to the boundary between the LWG and the LGM units.

Figure 6.

OBS 19 shear wave splitting observations, shallow anisotropic layer. (a) R and (b) T components. (c) Fast and slow components generated by rotation of X and Y using the symmetry plane orientations observed in (b). Left-hand panel is the rotated X component and right-hand panel is the Y component, using the symmetry plane angle which falls in the range 0–90°. (d) R and (e) T components following application of a layer stripping correction using the parameters derived in (b)–(c). (f) Location of instrument (inverted red triangle) and shots used to construct azimuthal stacks shown in (a)–(e) (blue lines) relative to OBS array (inverted grey triangles) and all GI shot lines (black lines) (see Fig. 2c). OBS data are filtered using a zero-phase low-cut (20–30 Hz) to remove a strong low frequency component at ∼9 Hz. Black arrows above and below R and T component panels show the mean orientation of the X and Y horizontal geophone components, and on fast and slow panels demarcate the edges of the 180°-wide polarity pairs (see Figs 1e and f). Dashed horizontal blue lines on R components provide indication of ‘flatness’ of arrivals. Dashed vertical blue lines on T components indicate locations of polarity reversals/azimuth nulls. Dashed and dotted light blue lines on fast and slow components show polarity pairs used for traveltime delay analysis. Numbered circles indicate locations of features discussed in the text. m indicates P–S time of direct-wave first multiple. Approximate stratigraphic correlations are labelled using the abbreviations given in Section 2.2.

Using the symmetry plane orientations of 30/120° to separate into fast and slow components and matching like-polarity pairs (Fig. 6c, events 3 and 4, comparable to Figs 1e–g, event 4) a time delay between the two components of ∼4–5 ms can be measured. Assuming a uniform accumulation of delay through all of the stratigraphy above the conversion point, this corresponds to 2–2.5 per cent anisotropy. The fast direction is oriented at 30°. Application of a layer stripping correction using an intermediate time delay of 4.5 ms results in improved flattening of the R arrivals (Fig. 6d, event 4, comparable to Fig. 1i event 6) and removal of the characteristic 4 × 90° polarity patterns from the T component (Fig. 6e, event 5, comparable to Fig. 1k, event 7).

Following the application of this shallow correction, analysis of deeper sources of anisotropy is performed. Inspection of the T component (Fig. 7b) shows that the amplitude nulls at 30/120° appear continuous down the record section. Several prominent azimuthally continuous arrivals are present on the R component at 0.57–0.60 s, 0.63–0.70 s and 0.84–0.96 s, with a further number of lower amplitude arrivals also present (Fig. 7a). These strong arrivals all occur at times consistent with P–S conversion points in the well-stratified AbG. Considering the first pair of high-amplitude R arrivals identified above (event 1, between 0.57 and 0.70 s), the corresponding time interval of the T component displays amplitude nulls with polarity reversals with an orientation of 20–30/110–120° (Fig. 7b, event 2), similar to that identified previously. Retaining the symmetry plane orientation of 30° and separating into fast and slow components (Fig. 7c) we find that arrival event 3 could in theory be paired with either event 4 or event 5. The corresponding time shifts in each case are 8 and 20 ms downwards and upwards, respectively. Applying the downward shift to match events 3 and 4 would result in a change of 90° in the orientation of the fast and slow directions. While such a change is possible under some circumstances, for example where the fractures generating anisotropy are pressurized with fluid (e.g. Crampin & Peacock 2005), we do not favour this interpretation. Instead we correlate events 3 and 5, with the 20 ms arrival time difference corresponding to a percentage anisotropy of ∼5 per cent if accumulated uniformly. For this solution there is no change in orientation of the fast and slow directions.

Figure 7.

OBS 19 shear wave splitting observations, deeper anisotropic layers. (a) R and (b) T components. (c) Fast and slow components generated by rotation of X and Y using the symmetry plane orientations observed in (b). Left hand panel is the rotated X component and right hand is the Y component, using the symmetry plane angle which falls in the range 0–90°. (d) R and (e) T components following application of a layer stripping correction using the parameters derived in (b) and (c). (f) Fast and slow components generated by rotation of X and Y using the symmetry plane orientations observed in (e). OBS data are filtered using a zero-phase bandpass filter (20–30–60–80 Hz). Annotation is as for Fig. 6.

Application of this value as a further layer stripping correction does appear to reduce the energy on the T component (Fig. 7e, event 6) at the corresponding time interval. Following this correction amplitude nulls can still be observed in the region 0.82–0.96 s (event 8, with corresponding energy on the R component, Fig. 7d, event 7) with the orientation of the symmetry planes of ∼30–120° matching those of the event above. Separating the fast and slow components and matching polarity-pairs (Fig. 7f, event 9) shows that a further time delay of 8 ms is present, corresponding to a percentage anisotropy of ∼4 per cent, with 30° continuing to be the fast direction.

6.1.2 OBS 20

On the T component of OBS 20, between 0.12 and 0.16 s, amplitude nulls and polarity reversals are observed at an orientation of ∼20/110° (Fig. 8b, event 2). Azimuthally continuous R component arrivals are also observed at this interval (Fig. 8a, event 1) and around 0.06–0.08 s (event 3). However, this shallower R component arrival does not display a corresponding 4 × 90° pattern on the T component with there instead being only a single polarity reversal at ∼190°.

Figure 8.

OBS 20 shear wave splitting observations, shallow anisotropic layer. (a) Rand (b) T components. (c) Fast and slow components generated by rotation of X and Y using the symmetry plane orientations observed in (b), ordered as in Fig. 6. OBS data are filtered using a zero-phase bandpass filter (20–30–60–80 Hz). (d) Location of instrument (inverted red triangle) and shots used to construct azimuthal stacks shown in (a)–(c) (blue lines) relative to OBS array (inverted grey triangles) and all GI shot lines (black lines). Annotation is as for Fig. 6.

Applying a rotation to separate into fast and slow oriented components and matching like-pairs (Fig. 8c, event 4), the observable delay between the two arrivals appears to be negligible at the level of the ∼1 ms picking precision that is possible, though it must still be sufficient to generate the SWS pattern of nulls on the T component. This small time delay is supported by the observation that the corresponding R arrivals appear to be essentially flat. Therefore, the percentage anisotropy associated with this arrival is low, but must be non-zero. For example, if a 2 ms delay was present the percentage anisotropy would be approximately 1.5 per cent, so we can infer the magnitude is less than this.

Looking instead at the first 1 s of the OBS records to attempt to identify evidence for deeper anisotropy, we see evidence for P–S converted arrivals on the R component (Fig. 9a) in the form of bright, azimuthally continuous arrivals at ∼0.30–0.34 (event 1), 0.61, 0.66–0.70 (event 3) and 0.86–0.96 s (event 4). There is a further strong azimuthally continuous arrival at ∼0.5 s (event 2), but the upwards curvature that parts of this arrival display suggest that it may either represent a P-wave arrival or a direct wave multiple. The normal moveout applied in this region should exclude the former attribution, and the arrival time does not fit well for the latter. This arrival is excluded from further analysis.

Figure 9.

OBS 20 shear wave splitting observations, deeper anisotropic layers. (a) R and (b) T components. Dashed grey lines in (b) are the symmetry plane orientation derived in Fig. 8. (c)–(d) as (a)–(b), focusing on arrivals in interval 0.6–1.0 s. (e) Fast and slow components generated by rotation of X and Y using the symmetry plane orientations observed in (d), ordered as in Fig. 6. (f) R and (g) T components following application of layer stripping using parameters from matching events 9 and 10 in (e). (h) Fast and slow components from (d) and (e) for the interval containing events 7 and 8. See discussion in text for further detail. OBS data are filtered using a zero-phase bandpass filter (20–30–60–80 Hz). Annotation is as for Fig. 6.

Few of these R component candidate P–S arrivals display a corresponding T component HTI SWS pattern (Fig. 9b). There are regions where parts of the characteristic patterns may be observed, and where either one 90° quadrant does not show clearly or the spacing does not appear to be exactly 90°. The clearest apparent T component arrivals appear to be those at 0.69–0.75 s, ∼30 ms after the corresponding R arrival (event 3). For the strong relatively flat R arrival event in box 5 in Fig. 9(c) there is no convincing evidence of a corresponding T component splitting pattern, which is instead observed at 0.7–0.75 s (Fig. 9d, event 6). Comparing the corresponding time interval on the R component shows that the arrival there is azimuthally discontinuous, which may suggest a relatively large time shift. At times >0.86 s we once again observe bright energy, and the R arrivals are undulating but remain azimuthally continuous (Fig. 9c, event 7). We also observe candidate T component HTI SWS patterns corresponding to these arrivals.

We apply a trial rotation of 70°, which is a compromise between the 60° and 75–80° angles that can be measured, to separate the fast and slow arrivals of the event at 0.69–0.75 s. This results in the fast and slow components shown in Fig. 9e). Matching of like polarity-pairs indicates a match between events 9 and 10 and gives a time shift of ∼10–12 ms and an orientation of the fast direction of ∼70°. Application of this as a layer stripping correction removes the transverse energy from box 6 (Fig. 9g) and transfers it upwards to box 5. The resulting R component arrival in box 6 is also now close to flat (Fig. 9f). However, the post-layer stripping R component arrivals in box 5 and at events 7 and 8 (Fig. 9f) now display significant deterioration in their apparent flatness. Furthermore, the locations of the peaks and troughs of these undulating reflectors implies that the fast and slow directions are 90° different to those for event 6, that is they are transposed. As discussed in the previous section, while it is possible to generate this transposition through, for example, fluid pressure, it is highly unlikely that this process is acting in this location, where there is no observed gas, and over the narrow stratigraphic interval required by the proximity of these observations in time. It is, therefore, unlikely that event 6 is a primary P–S arrival displaying HTI-type SWS.

If we ignore the events at 5 and 6 and analyse those at 0.86–0.96 s (Figs 9d and e), we can see that both the strong R component arrivals at 0.86–0.89 s (event 7) and 0.92–0.95 s (event 8) display a corresponding T component with polarity reversals spaced at ∼90°, and an orientation of the symmetry planes of ∼60–70° (as above for event 6). Separation into fast and slow components (Fig. 9h) and matching of the polarities marked events 11 or 12 both give a traveltime delay of ∼4 ms and a fast direction of 70°. Application of this correction does a good job at removing the small observable cos2θ traveltime variation observed on events 7 and 8 (Fig. 9d). Due to the relatively complex nature of the SWS observations, and the fact that they are an integrated measurement over the ray path, it is not possible to quantify the strength of the anisotropy, as it is difficult to identify the interval over which the traveltime delay accumulates. However, it is likely that the value is small, on the order of a few per cent.

6.2 Scanner Pockmark, OBS 1

OBS 1 is located within West Scanner Pockmark (Figs 2c and 10f). Observations of shallow anisotropy on this instrument were presented by Bayrakci et al. (2021). Following the adjusted methodology applied in this paper to target deeper arrivals than was possible in that previous study, we also observe symmetry planes oriented at ∼70°/160° (Figs 10a, b events 1 and 2). Using these values for separation into fast and slow components and matching like-polarity pairs (Fig. 10c, events 3 and 4) results in a range of measured delay values from 0 ms (i.e. no measurable delay) up to an upper bound of 2 ms and a corresponding fast direction of 160°. Bayrakci et al. (2021) did not find a reliable value for the traveltime delay and concluded that the small magnitude of this may arise due to the very short path length that does not allow the accumulation of a delay between the fast and slow S-waves above the level of picking error. Our observation of an upper bound of 2 ms delay is consistent with this inference. Based on the P–S time after the direct arrival, if this delay is present then it would correspond to ∼2.5 per cent anisotropy. By applying the upper bound of 2 ms as a layer stripping correction we see that there is a removal of the characteristic T component pattern (Fig. 10e, event 6). Due to the small magnitude of the time delay it is difficult to resolve any change in the R component pattern (Fig. 10d, event 5). Overall, our observations of shallow anisotropy support those made previously by Bayrakci et al. (2021).

Figure 10.

OBS 1 shear wave splitting observations, shallow anisotropic layer. (a) R and (b) T components. (c) Fast and slow components generated by rotation of X and Y using the symmetry plane orientations observed in (b), ordered as in Fig. 6. (d) R and (e) T components following application of a layer stripping correction using the parameters derived in (b) and (c). OBS data are filtered using a zero-phase low-cut filter (20–30 Hz). (f) Location of instrument (inverted red triangle) and shots used to construct azimuthal stacks shown in (a)–(e) (blue lines) relative to OBS array (inverted grey triangles) and all GI shot lines (black lines). Dashed orange lines indicate extent of seabed expression of West and East Scanner Pockmarks. Annotation is as for Fig. 6.

Following the analysis of the shallow anisotropy and application of the upper bound 2 ms layer stripping correction, we look for evidence of SWS in deeper layers. Clear azimuthally continuous arrivals at ∼0.7 and 0.81 s (Figs 11a and c) indicate P–S conversions are occurring within the well-stratified AbG. Between the events at ∼0.1 and ∼0.7 s there are no coherent events except for the direct arrival multiple at 0.18–0.22 s. Instead there is a chaotic and disorganised pattern of arrivals, which is mirrored in the T component (Fig. 11b). Stratigraphic correlation indicates that this region spans the CP and up to 85 m of the combined thickness of the LB and AbG. The uncertainty on this latter value results from a lack of information about the degree of gas saturation in the chimney, and, hence, the extent of velocity suppression there. As gas presence would lower velocities in this region, the 0.7 and 0.81 s events may correspond to shallower conversion depths, though these would still be expected to be within the AbG.

Figure 11.

OBS 1 shear wave splitting observations, deeper anisotropic layers. (a) R and (b) T components. (c) and (d) as (a) and (b), focusing on arrivals in interval 0.55–1.1 s. Dashed grey lines are the symmetry plane orientation derived in Fig. 10. (e) Fast and slow components generated by rotation of X and Y using the symmetry plane orientations observed in (d), ordered as in Fig. 6. OBS data are filtered using a zero-phase bandpass filter (20–30–60–80 Hz). Annotation is as for Fig. 6.

The R component arrival at 0.7 s does not show any corresponding clear HTI-type patterns on the T component (Figs 11c and d, event 2). However, the arrival at 0.81 s (event 3) shows evidence for polarity reversals suggesting symmetry plane angles of ∼15/105°, although the pattern is unclear or missing in the 180–270° quadrant and the azimuths may be close to that of the shallow anisotropy (Fig. 11d, grey dashed lines). Using 15/105° as the orientations of the symmetry planes results in good separation into fast and slow components (Fig. 11e). Matching of like-polarity pairs gives the smallest time-shift solution of 2–4 ms (event 4), with 15° as the fast orientation. Due to the low frequencies and irregularity of the separated phases increasing the picking error in this region, and the uncertainty over the symmetry plane orientations, this value likely represents an upper bound. If this delay is assumed to accumulate between ∼0.7 and ∼0.81 s it would correspond to ∼2–4 per cent anisotropy, which is similar to that observed on the OBSs at the reference site (≤4–5 per cent).

6.3 Instruments around pockmark

6.3.1 OBS 8

OBS 8 is located immediately to the south of West Scanner Pockmark, ∼200 m from the pockmark's centre (Figs 2c and 12d). This instrument was also studied by Bayrakci et al. (2021). Utilizing the different approach to flattening P–S arrivals that is applied in this study, we make the same observations of anisotropy as the previous study. The symmetry planes are observed to be oriented at ∼70/160° (Fig. 12b, event 2). The corresponding R component arrival shows some azimuthal variation with time (Fig. 13a, event 1) that may be attributed to small residual variability of the direct arrival time, which is not fully corrected during all of the processing steps. Using these symmetry plane orientations, rotation into fast and slow components and like-polarity pair matching (Fig. 12c, events labelled 3) shows that the measurable time delays are no greater than the pick uncertainty. Bayrakci et al. (2021) attributed this observation to a lack of opportunity for the signals to accumulate a delay, given that the arrival is only 70–80 ms after the direct arrival. Given it should be possible to pick delays at this depth with an uncertainty of 1–2 ms we can infer that the anisotropy has a magnitude of <∼1.5–2 per cent.

Figure 12.

OBS 8 shear wave splitting observations, shallow anisotropic layer. (a) R and (b) T components. (c) Fast and slow components generated by rotation of X and Y using the symmetry plane orientations observed in (b), ordered as in Fig. 6. OBS data are filtered using a zero-phase low-cut filter (20–30 Hz). (d) Location of instrument (inverted red triangle) and shots used to construct azimuthal stacks shown in (a)–(c) (blue lines) relative to OBS array (inverted grey triangles) and all GI shot lines (black lines). Dashed orange lines indicate extent of seabed expression of West and East Scanner Pockmarks. Annotation is as for Fig. 6.

Figure 13.

OBS 8 shear wave splitting observations, deeper anisotropic layers. (a) R and (b) T components. (c) Fast and slow components generated by rotation of X and Y using the symmetry plane orientations observed in (b) for the interval 0.43–0.52 s, corresponding to location of event 2, and ordered as in Fig. 6. (d) Fast and slow component separations for interval 0.81–0.92 s, corresponding to event 4, following application of a 3 ms time shift to correct for anisotropy associated with event 2 and shown in (c). OBS data are filtered using a zero-phase bandpass filter (20–30–60–80 Hz). Annotation is as for Fig. 6.

At greater depths there are several clear sets of azimuthally continuous R component arrivals (Fig. 13a), with the highest amplitudes at 0.45–0.50 (event 2) and 0.82–1.0 s (event 4). Stratigraphic correlation places these events in the AbG. Weaker candidate arrivals are also visible at ∼0.28 s (event 1) and ∼0.66 s (event 3). At the corresponding time intervals on the T component (Fig. 13b), characteristic patterns associated with HTI-type SWS are observed at 0.25–0.30, 0.45–0.50 (although the 0–90° azimuth quadrant pattern is less clear here) and 0.82–1.0 s. In all of these cases the orientation of the symmetry planes appears to be ∼45/135°, which is close to the regional stress direction of ∼54/144° (Evans & Brereton 1990). We exclude events 1 and 3 from further analysis because the former may be affected by low-frequency interference from the direct wave immediately above, and the latter does not have a corresponding characteristic T component SWS pattern.

For the remaining events 2 and 4, we apply a fast and slow component separation using the symmetry plane angles of 45/135° (Figs 13c and d), based on the observations from the T component. For event 2 the separation angle appears to be suboptimal, as there is some smearing between the 180° separated arrival pairs. Using this value we measure a time delay of 3–5 ms (Fig. 13c, events 5 and 6) and a fast direction of 45°. If it is assumed that traveltime delay is accumulated uniformly below the arrivals at 0.06–0.09 s the resulting percentage anisotropy is 0.8–1.4 per cent, which is consistent with that suggested for the upper layer. Event 2 on the R component shows a lateral break in continuity between the 180–270° and 270–360° quadrants, while there is a lack of a clear break in the T component arrivals in the 0–90° quadrant. These parts of the shot-receiver azimuth range cover the northern side of the instrument. The location of OBS 8, ∼200 m south of the pockmark centre (Fig. 12d), and the short offset range of the shots used in constructing this azimuth stack means that seismic energy may be disrupted by the nearby gas and/or chimney structure. Based on the shot and receiver positions (Fig. 12d), and since Vp >> Vs, we expect any disruptive effects due to the presence of the chimney and gas to occur during the downgoing P-wave leg.

Following the application of a 3 ms time-shift to correct for the anisotropy observed in the interval 0.45–0.5 s, subsequent analysis of the corresponding separated components at ∼0.86 s (Fig. 13d) using the same orientation of the symmetry planes shows that there is no robustly measurable delay above the picking error for the frequencies used (event 7). This again indicates a very small degree of anisotropy. For example, if the magnitude of anisotropy were the same as in the interval 0.1–0.45 s found above, then the expected range of delay times would be ∼3–5.5 ms, which should be observable at this level.

6.3.2 OBS 12

OBS 12 is located ∼400 m ESE of the West Scanner Pockmark and ∼250 m SE of the smaller East Scanner Pockmark. There are two different sets of profiles which can be used to generate an approximately equal offset square of shots around the instrument (Fig. 14f). Comparison with the radial component semblance (Fig. 4b) shows that for shallow arrivals the shorter offset shots are preferable since there is greater localization of semblance peaks, indicating a greater coherence of the seismic energy in this region. In contrast, for the deeper arrivals we observe a degradation in the velocity localization of the semblance peaks, indicating that for these regions we need to use arrivals at longer offsets.

Figure 14.

OBS 12 shear wave splitting observations, shallow anisotropic layer. Shots used in constructing seismic sections shown are short offsets, shown by the inner blue rectangle in (f). (a) R and (b) T components. (c) Fast and slow components generated by rotation of X and Y using the symmetry plane orientations observed in (b), ordered as in Fig. 6. (d) R and (e) T components following application of a layer stripping correction using the parameters derived in (b) and (c). OBS data are filtered using a zero-phase low-cut filter (20–30 Hz). (f) Location of instrument (inverted red triangle) and shots used to construct azimuthal stacks shown in (a)–(e) (blue and purple lines) relative to OBS array (inverted grey triangles) and all GI shot lines (black lines). Dashed orange lines indicate extent of seabed expression of West and East Scanner Pockmarks. Annotation is as for Fig. 6.

Using the short offset shots (<∼150 m, inner blue square on Fig. 14f), we observe azimuthally discontinuous arrivals on the R component in the time interval 0.14–0.18 (Fig. 14a, event 1). In the corresponding time interval on the T component, azimuth nulls can be observed with an orientation of 60–70°/150–160° (Fig. 14b, event 2). Using 70° as the angle of the symmetry plane to separate into fast and slow components, a delay of 3–5 ms is observed between like-polarity pairs (Fig. 14c, event 3). This corresponds to a percentage anisotropy of ∼2.1–3.5 per cent and a fast direction of 70°. Application of a mean value of 4 ms as a layer stripping correction removes the lateral discontinuities in the R component arrivals (Fig. 14d, event 4). Stratigraphic correlation places the conversion from P to S towards the middle of the Lower Witch Ground Fm., which sub-bottom profiler imaging shows to be well layered (Fig. 3b).

In the deeper layers, we can observe the effects in the R and T component azimuth stacks of using longer (∼300–450 m; Figs 15a and b; outer purple profiles on Fig. 14f) and shorter (<∼150 m; Figs 15c and d; inner blue profiles on Fig. 14f) shot offsets. For both sets of offsets there are flattened arrivals in the interval 0.3–0.35 s (Fig. 15a event 1 and Fig. 15c event 3). For the shorter offsets (event 3) these have a high amplitude and appear to merge upwards into the latter parts of the direct arrival multiple at ∼0.2 s. This can be seen in the merging of peaks between the P–S semblance trend and the direct arrival multiple in the short offset semblance (Fig. 5b, 50–200 m). It is possible, therefore, that some of the energy in this region is not related to P–S conversions. Furthermore, there does not appear to be any clear evidence for HTI-type SWS patterns in the corresponding region of the T component (Fig. 15d, event 4). Looking instead at the longer offsets (200–400 m), the R component semblance shows some weaker but isolated peaks in this region (Fig. 4b). The R component azimuthal stack (Fig. 15a) shows evidence for an azimuthally continuous arrival (event 1) which contains two small breaks at ∼180° intervals. The corresponding T component shows some evidence for polarity reversals, particularity in the range 180–360°, and which give a direction of anisotropy of ∼30/120° (event 2). However, no clear fast/slow phase separation can be conducted and so anisotropy cannot be confirmed or quantified for these arrivals.

Figure 15.

OBS 12 shear wave splitting observations, deeper anisotropic layers. (a) R and (b) T components, longer offset shots (outer purple profiles, Fig. 14f). (c) R and (d) T components, shorter offset shots (inner blue profiles, Fig. 14f). (a)–(d) OBS data are filtered using a zero-phase bandpass filter (20–30–80–100 Hz). (e) R and (f) T, longer offset shots, filtered using a zero-phase bandpass filter (20–30–60–80 Hz). (g) Fast and slow components generated by rotation of X and Y using symmetry plane orientations derived using (e) and (f), ordered as in Fig. 6. Annotation is as for Fig. 6.

At greater P–S times, corresponding to conversions within the AbG, packages of azimuthally continuous, near-flat arrivals are observed on the R component at 0.84–0.92 and 0.95–1.05 s (Fig. 15e, events 5 and 6), for the longer-offset (300–450 m) shot profiles. The upper of these events does not show any HTI-type SWS patterns at the corresponding time on the T component, but the lower event does (Fig. 15f, event 7), with polarity reversals observed with an azimuth of ∼50/140°. Application of the rotation to fast and slow components reveals the presence of separate phases (Fig. 15g, event 8), with a time delay on the order of ∼2–4 ms, and giving a fast direction of 140°. This is consistent with the already flat appearance of the arrival on the R component. In addition, given the flat appearance of the earlier 0.84–0.92 s radial arrivals, we infer that the observed anisotropy accumulates only in the interval between the two arrival packages. Given the ∼0.1 s time-span of each package, the corresponding magnitude of anisotropy is uncertain, but if the top of both is used to estimate this interval we would interpret a magnitude of anisotropy of ∼2–4 per cent.

6.4 Summary of observations

Observations of seismic anisotropy have been made on five OBSs in and around the Scanner Pockmark Complex. One instrument was located within the pockmark itself, with another located on the southern edge of the pockmark rim, and three further OBSs at a distance of ∼400 m to 1.5 km from the pockmark centre (Fig. 2c). Fig. 16(a) summarizes the observations made, highlighting the anisotropic symmetry plane orientations, the fast directions, and the magnitude of anisotropy, where they are possible to derive. There is variation in both the spatial and depth extent of anisotropy in and around Scanner Pockmark, which we relate to the different subsurface structures present in and around the pockmark (Fig. 16b).

Figure 16.

Summary of SWS observations in and around Scanner Pockmark Complex. (a) Summary of anisotropic symmetry plane orientations, with fast direction orientation where available and percentage anisotropy labelled for stratigraphic units in the Scanner Pockmark study area. Values are italicized where considered less robust. Vertical bars indicate approximate extent of anisotropic regions, assuming uniform accumulation of traveltime delay. (b) Schematic diagram summarizing relationships between subsurface structures within (OBS 1) and on flank (OBS 8) of Scanner Pockmark and at away from pockmark and at reference site (OBSs 12, 19 and 20). Stratigraphic annotations as for Fig. 3.

Observations of anisotropy made in this study are limited to those where HTI-type SWS patterns are visible, as the layer stripping procedure to successively correct shallow anisotropic layers is only applicable in this case. In calculating the percentage anisotropy for various layers and observations, we have assumed that the time delay between the fast and slow arriving shear components is accumulated uniformly. However, this may not necessarily be the case, because there may be unresolved layer boundaries. Therefore, our estimates for percentage anisotropy, particularly for the deeper subsurface layers, represent lower bounds.

There is uncertainty in the interpretation of which layers generate the anisotropy which may affect the interpretations of sources of anisotropy discussed below. The approach used to develop the subsurface velocity structure in this study is highly sensitive to very thin layers present in the stratigraphy. Uncertainties in Vp and Vs in the very shallowest subsurface may result in uncertainties of 10 s of milliseconds in the P–S time of the converting horizon, particularly for the shallowest layers. However, this uncertainty is unlikely to impact significantly on the interpretations of shallow anisotropy, as these observations still primarily correlate to structures at the base of the Upper Witch Ground or within the Lower Witch Ground. The observations of deeper anisotropy are predominantly found at a P–S time of ∼0.7 s below the seabed, which, even accounting for the uncertainty in the WG unit velocities, places the conversions from P to S within the upper parts of the AbG. This conversion depth is consistent with the well-layered seismic appearance of this formation (Figs 3d–f).

7 DISCUSSION

7.1 Regional anisotropy

The regional, or ‘background’ pattern of anisotropy is most evident away from the pockmark in the Aberdeen Ground Fm., which displays evidence for seismic anisotropy throughout our study area. This unit generates clear conversions from P- to S-waves, observed as azimuthally continuous arrivals on the radial components (e.g. Figs 7 and 9), due to its internally layered seismic structure (Figs 3d–f). The observations of anisotropy on OBSs 12, 19 and 20 in this layer are not uniform, with the strength varying between ∼1 and 4 per cent and a ∼40° range of the generally NE–SW fast direction (Fig. 16a).

Strong P- to S-wave conversion also appears to occur between the Upper and Lower Witch Ground Fms., Lower Witch Ground and LGM Fms., and at the top of the Aberdeen Ground Fm.. Away from the pockmark conversions are also occasionally observed between the LGM and Coal Pit Fms.. However, the boundaries between other layers are not always strong, so we might not expect these to produce clear P–S converted arrivals. For the UWG–LWG and in some cases the LWG–CP boundaries, the shallow depth means that there may be insufficient S-wave traveltime to generate a resolvable delay, particularly if the magnitude of anisotropy is on the order of the 1–2 per cent indicated by this study (Fig. 16a). For the intermediate depth range in the study area, comprising the LGM and CP Fms., the identification of potential anisotropy is limited by a lack of clear reflection boundaries in the seismic reflection data (Figs 3d–f) and the survey geometry which results in a strong direct arrival multiple overprinting the signals in the corresponding time interval.

Direct comparisons between our observations and previous studies of anisotropy are challenging as often very different targets are being investigated. However, previous studies of seismic anisotropy using SWS of active source seismic data have targeted a range of sedimentary and tectonic environments, which can be used to compare against the magnitudes of anisotropy found in this study. In a study located on the continental slope west of Svalbard, Haacke et al. (2009) find anisotropy of 1–2 per cent in a region of shallow (∼500 m thickness) sediments where accumulations of gas hydrate and free gas had previously been detected. Similar magnitudes of ∼2 per cent were also found within the uppermost 35 m of Quaternary deposits in the central New Madrid Seismic Zone (Harris 1996). Larger values of between 4 and 14 per cent were obtained across a range of depths within two oil fields located within the San Joaquin Basin, California (Winterstein & Meadows 1991). Overall, therefore, our values of ∼1–4 per cent anisotropy observed around the Scanner Pockmark region are consistent with the magnitudes of anisotropy observed using similar approaches in other locations and geological settings.

7.2 Implications for the understanding of chimney structures

7.2.1 Structure within the chimney

The structure of the physical fluid flow chimney will significantly impact our ability to resolve its properties using the techniques applied in this study. Several different geometries for the fluid flow network permitting vertical migration of gas through the physical chimney could be considered (e.g. Bull et al. 2018; Roche et al. 2021). If fluid flow occurs along vertically or subvertically oriented linear features then it would be expected to see evidence for HTI-type anisotropy. These linear features may correspond to pre-existing preferred orientations or regional stresses. Alternatively, the physical chimney could be comprised of a network of fractures that may not be vertically aligned, and might instead be dominantly concentric or radial. Over a spatially restricted region it would be possible to consider and model each of these idealized structures as an HTI anisotropic system, but in either case the fracture orientation will vary across the actual feature. Therefore, depending on the spatial scale of the target, experiment and data used, the effects of small-scale apparent HTI could be cancelled out. Finally, the sedimentary material within the chimney structure may itself be different to that which surrounds it. This may be the case if, for example, an initial pockmark formation blow-out event was followed by infilling with material of different porosity and permeability. In this case, we would not expect to see fracture-driven anisotropy, since the fluid migration is not focused by fractures. Alternatively, precipitation of methane-derived authigenic carbonates inside fractures may reduce the ability of the fractures to transmit fluids (e.g. Hovland et al. 1987).

The lack of observed P–S arrivals at 0.1–0.7 s beneath the pockmark (OBS 1, Figs 11 and 16a) is likely to arise both as a result of overprinting by the direct arrival multiple between 0.16 and 0.23 s, and as a result of an acoustic blanking effect due to the accumulation of gas in a shallow reservoir immediately below the pockmark (Bayrakci et al. 2021; Callow et al. 2021) and in the LB at the top of the seismic chimney (Figs 3d and e). However, it is also possible that physical chimney formation processes (erosive fluidization and blow-out) alter the stratigraphy within it such that converting boundaries are not present in this region. The re-appearance of flattened R component arrivals at ∼0.7 and 0.81 s corresponds to a depth below the top of the LB gas accumulation of up to ∼85 m, with this being an upper bound since velocities in this region may be reduced by the presence of gas (Schramm et al. 2021). These conversions, therefore, originate from stratigraphic horizons within the AbG. The P–S events at 0.81 s indicate an upper bound of ∼2–4 per cent anisotropy and a fast direction of ∼15°, close to N–S, if it is assumed that all the observed traveltime delay is accumulated between ∼0.7 and 0.81 s. This would be consistent with the degree of anisotropy observed in the AbG at the reference site (e.g. Fig. 7; Section 6.1.1). Alternatively, if the traveltime delay is accumulated over a greater stratigraphic thickness, obscured by the presence of the gas shadow region, then a smaller magnitude of anisotropy would be inferred.

Using a representative value for anisotropy based on observations outside the pockmark of 2–4 per cent (Fig. 16a) we can estimate the delay time which would be predicted to accumulate between 0.1 and the R component arrival at 0.7 s (Fig. 11c, event 2), that is through the region shown by Fig. 11(a), label 1. If this background level of anisotropy were present here, then over this interval we would expect to accumulate up to ∼10–20 ms of traveltime delay. A delay of this magnitude should be readily observable in the data, so we consider its presence to be unlikely. Thus, we suggest that the degree of HTI-type anisotropy within the fluid migration system within the pockmark is lower than that outside the chimney.

The ∼15° orientation of the fast SWS component (Fig. 16a) is similar to the anisotropy in the LGM unit above the gas (160/340°, Figs 10 and 16a). This is also similar to the direction of the long axis of both the West Scanner Pockmark (e.g. Fig. 2c) and the small pockmarks throughout the region (Böttner et al. 2019). This similarity suggests that the causative mechanisms generating the observed anisotropy and the formation of the pockmark are likely to be related. Our observations suggest that the initial pockmark blow-out results in an overprinting of the original ‘background’ anisotropy, with the post-blow-out chimney structure displaying ∼N–S oriented HTI-type anisotropy. This anisotropy is most likely to be relatively small in magnitude, up to ∼2 per cent.

7.2.2 Structure outside chimney

OBS 8, located immediately to the south of the West Scanner Pockmark, displays evidence for HTI-type anisotropy over three different depth intervals. The shallowest anisotropy is interpreted by Bayrakci et al. (2021) as relating to very shallow (<5 m depth) vertically aligned gas conduits, which facilitate the ebullition of methane at the seafloor. These features were interpreted as resulting either from the opening of fractures perpendicular to the regional minimum horizontal stress or due to the local stress gradient driven by overpressure arising from the gas accumulation in the LB.

The two deeper observations of anisotropy beneath OBS 8, at 0.45–0.50 s and >0.82 s, display a 45/135° orientation of the symmetry planes (Fig. 16a). Stratigraphic correlation places the 0.45–0.5 s event close to the top of the seismic chimney structure. Thus, the traveltime delay observed here is accumulated in the LGM and CP. Fig. 3(d) shows that there is a bright reflection at the top of the seismic chimney, immediately adjacent to the north, indicating the extent of the LB gas accumulation. Below this we can see that the seismic chimney structure does not extend beneath OBS 8. This may explain why many P–S conversion events are present at times >0.45 s, as this layer is the well-stratified and reflective AbG. The anisotropy observed in the 0.82 s event is inferred to be accumulated within the AbG and/or part of the LB. The percentage anisotropy associated with both depth intervals is small, at ∼1 per cent, with an observed fast direction of 45° (Fig. 16a). This direction is several 10 s of degrees different from the directions of anisotropy observed within the seismic chimney, and is more oriented towards the direction of the background anisotropy, which is dominantly NE–SW oriented (Fig. 16a). This observation suggests that this location is not directly influenced by the formation of the physical fluid-flow system. However, the lower magnitude of the anisotropy observed on this instrument (∼1 per cent; Fig. 16a) compared to OBSs 12, 19 and 20 which are at a greater distance from the pockmark suggests that there is some remaining influence of the physical gas migration system at distances of ∼100–150 m from the pockmark centre (Fig. 16b).

7.3 Origin of anisotropy

A range of potential causative mechanisms could explain the observed anisotropy in and around Scanner Pockmark (Bayrakci et al. 2021; Callow et al. 2021). For example, we might expect tensional fractures to form at ∼140–150°, perpendicular to the regional minimum horizontal stress direction, σ3, of ∼54° (Fig. 16b; Evans & Brereton 1990). Such a fracture system is not consistent with the observed anisotropy at the SPC and reference site which has the fast direction (where this can be determined) oriented in an overall NE–SW direction, that is line with the maximum stress (Fig. 16). It is possible for linearly oriented fractures to display a fast direction of anisotropy perpendicular to the fracture direction, for example when the fractures are pressurized with fluid (e.g. Crampin & Peacock 2005) as could be the case within the physical chimney. However, this interpretation would not explain the observations of anisotropy on OBSs 19 and 20 where no gas migration system is observed to be present (Fig. 3f), and, hence, we do not favour it.

Analysis of 3-D seismic reflection data reveals extensive linear features across the study area (Bayrakci et al. 2021; Callow et al. 2021). Small linear features at the base of the WG, interpreted to represent iceberg ploughmarks, have a dominantly NE-to-E orientation in the range ∼50–80°. Larger features interpreted as mega-scale glacial lineations (MSGLs) are oriented in a more NW–SE direction, of ∼150–160° (Fig. 16b). Both these features are also visible in 2-D sub-bottom profiler data, with the ploughmarks having a thickness of only a few ms TWTT and the MSGLs having a thickness up to ∼15 ms TWTT, present within most of the thickness of the LGM unit (Callow et al. 2021). Within the LB, the gas accumulation shows two principal trends. Gas appears to be focused at orientations of ∼60° and ∼320° (or 140°), extending from the centre of the pockmark towards the north. Scanner Pockmark also lies on the northern edge of a glacial tunnel valley which forms part of the LB units, eroding into the underlying AbG (see Fig. 3e). This feature is oriented in an approximate NW–SE direction (Callow et al. 2021), with the edge of the tunnel valley passing close to OBS 8.

The consistency between the orientation of one of the directions of gas spreading away from the centre of the pockmark in the LB with the approximately NE–SW oriented fast directions on OBSs 12, 19 and 20 at corresponding stratigraphic intervals suggests there may be a correlation between the causative mechanisms. For example, this may indicate a trend of higher permeability in the unit beneath the seal, even though there is no gas beneath those instruments. This direction is also consistent with the orientation of the iceberg ploughmarks, which imply a flow direction governed by bottom water currents and/or tides. While the effect of the iceberg ploughmarks themselves on anisotropy is likely to be small, due to their limited vertical extent, they may help to indicate an overall control by the same larger-scale processes which generate greater anisotropy at depth. Therefore, these observations suggest that the ‘background’ anisotropy in the sediments results from syn-depositional and post-depositional glaciomarine processes. This interpretation may also explain why we observe a degree of variability in the orientation and strength of the background anisotropy, particularly in the Aberdeen Ground Fm..

With respect to the anisotropy beneath the pockmark, the approximate N–S orientation of the fast direction is consistent with the trend of the long axis of the West and East Scanner Pockmarks (Fig. 2c). This suggests a potentially linked causative mechanism between the two observations, and, thus, that the gas migration at depth in the present day is governed by the same controlling factors as the initial pockmark and chimney formation, (Cevatoglu et al. 2015; Callow et al. 2021; Roche et al. 2021). However, the strength of anisotropy in the region immediately below the pockmark remains unclear.

Finally, the anisotropy observed at the southern flank of the pockmark, on OBS 8, provides some indication of a lateral influence of the fluid flow processes within the chimney. This observation may indicate the limits of the spatial extent of the migration system, or may result from local stress gradients arising from overpressure in a more spatially restricted migration system.

8 CONCLUSIONS

We observe evidence of SWS associated with vertically aligned linear features (HTI-type anisotropy) on several ocean bottom seismometers located within and adjacent to the Scanner Pockmark Complex, and at a nearby reference site. By determining the orientations of the anisotropic symmetry planes, the magnitude and fast direction of anisotropy and the stratigraphic interval of origin we find that:

• Seismic anisotropy is commonly observed in arrivals originating from P-to S-wave conversions, particularly from the shallowest Witch Ground Formation and the deeper, well-layered Aberdeen Ground Formation. Overall, the magnitude of background anisotropy in the region is small, on the order of ∼1–4 per cent, with the largest values being observed in the Aberdeen Ground Fm. away from the pockmark.

• The region beneath Scanner Pockmark containing the physical shallow gas migration system and the seismic chimney below the gas accumulation in the Ling Bank Fm. both show evidence for anisotropy, which is oriented approximately N–S, in line with the long-axis of the pockmark depression, suggesting this anisotropy is related to the pockmark formation process. The observation of the deeper anisotropy suggests the presence of a deeper fluid migration system that feeds the shallow reservoir at the base of the physical chimney.

• The gas migration system beneath the pockmark displays a relatively small degree of anisotropy throughout, which indicates that the gas migration system is not strongly governed by vertically aligned anisotropy. However, since parts of this region are obscured either by the direct arrival multiple or by the effect of the gas shadow, only an upper limit on anisotropy can be defined here.

• Observed changes in the direction and strength of anisotropy at the edge of the pockmark may indicate the limits of fracturing within the chimney and/or the effects of overpressure gradients, and mark the bounds of the fluid migration system.

• There is a correlation between the lateral spreading of gas at depth in the Ling Bank Fm. and the orientation of ‘background’ anisotropy in the region, which we suggest is related to syn- and post-depositional sedimentary processes involved in the formation of the Quaternary sediment succession.

• Overall, there are differences in anisotropy associated with the gas migration system beneath Scanner Pockmark relative to the surrounding region. Thus, we conclude that this technique represents a promising tool for evaluating chimneys as potential leakage pathways. However, to further examine gas migration systems both above and below gas accumulations, considerations must be made to optimize the survey geometry in order to maximize the imaging potential across all depth intervals, and to better determine the spatial extents of the system.

ACKNOWLEDGEMENTS

This work was funded by the Natural Environment Research Council (CHIMNEY; NERC Highlight Topic; NE/N016130/1 and NE/N015762/1) and the European Union's Horizon 2020 research and innovation programme under grant agreement No. 654462 (STEMM-CCS). We would like to thank all those involved in the planning and acquisition of data during RRS James Cook cruise JC152, including the officers, engineers and crews, the scientific parties, and all seagoing technicians and engineers. OBSs were supplied by the UK Ocean Bottom Instrumentation Facility (Minshull et al. 2005). We are also grateful for the support of Applied Acoustics Ltd during sparker data acquisition. We thank Ross Haacke (CGG) for advice on the design of the CHIMNEY seismic surveys and for providing codes for rotation of OBS data and Giuseppe Provenzano for work processing of the seismic reflection data. We thank Martha Savage, Sunny Singhroha and editor Jenny Collier for their helpful reviews. The final accepted version of this manuscript is available through ePrints Soton (eprints.soton.ac.uk).

DATA AVAILABILITY

All project data have been archived by the British Oceanographic Data Centre (BODC). OBS data used in this study are available at http://dx.doi.org/10.5285/0068C20ACA6C43E786E737C9257179A2

APPENDIX A: OBS COMPONENT TILT AND ROTATION ANALYSIS

A.1 INTRODUCTION

As discussed in Section 4.1.2, during data QC it was noted that the X geophone component on several instruments appeared to show a delay to the peak of the direct arrival when compared to the Y and Z components (Figs A1a-c and A2a-c), accompanied by a longer wavelet, indicating a lower frequency content (Figs A1f and A2f). The function of these observed behaviours on the data appeared to take the form of a non-zero-phase filter, which for simplicity we separated into two components, a zero-phase low-pass filter and a static delay. Hereafter we refer to this as the X geophone functional filter.

Figure A1.

OBS geophone alignment correction, OBS 1. Corrections are applied to reduce or remove the effect of the X geophone functional filter. (a) X, (b) Y, (c) Z components, flattened on the direct arrival and centred at 0 ms using the amplitude of the Z component. Offsets shown are 200–220 m (absolute), and traces are ordered by shot number. The first peaks of the X component are delayed relative to Y and Z. (d) Xa component (re-alignment case) generated by application of a bulk time-shift to X, and (e) Ya, which is identical to Y since correction is applied to X. Improvement in the alignment of peaks between the two horizonal components can be seen. (f) Frequency spectra of stacked X (red) and Y (blue) traces in 100–300 m offset range. Panel (g) Xa,f and (h) Ya,fre-alignment + filter horizontal geophone components.

Figure A2.

OBS geophone alignment correction, OBS 19. Plotted as for Fig. A1. Note there is little to no delay observed between the X, Y and Z components (a–c), and therefore additional re-alignment correction (d–e) has little effect. The frequency spectra (f) of the X and Y components are also very similar, in contrast to Fig. A1, and so the re-alignment + filter correction also shows virtually zero effect (g–h).

Figure A3.

Tilt and rotation analysis, OBS 1. Tilt analysis outputs for (a) initial, (b) re-alignment and (c) re-alignment + filter test cases. α and β are the X and Y tilt angles respectively, where values of 90° for both indicates these components are perfectly horizontal. White star indicates the location of the calculated misfit minimum, with value labelled. (d) Calculated trace-by-trace rotation angles for the initial, (a, black), re-alignment (b, red) and re-alignment + filter (c, blue) cases prior to application of tilt-correction. Panel (e) as (d), for post-tilt corrected components. Blue line shows the fit to the blue circles, equation of this curve (R_out) and r² fit value are labelled. (f) Output trace-by-trace radial:transverse power ratio for the pre-tilt corrected cases, colours as in (d). Panel (g) as (f), for post-tilt corrected components.

Figure A4.

Tilt and rotation analysis, OBS 19. Plotted as for Fig. A3. Note that in (a–c), the value of the standard deviation misfit minimum is smaller than for OBS 1 (Fig. A3), where the X geophone functional filter was stronger. In panels (d) and (e), note that the peak-to-peak amplitude variation of the rotation angle around the mean values is also decreased compared to OBS 1, and the three test cases show a smaller overall range. In panels (f)–(g), note how the R/T power for the initial test case is higher than for OBS 1.

The presence of the X geophone functional filter impacts processes which depend on the alignment and/or measurement of properties of the direct arrival. Two key steps in the data analysis process affected by this are (1) geophone package tilt analysis and correction and, (2) rotation to R-T components. To correct for the X geophone functional filter, we apply a static shift to the X geophone component, based on the average delay measured for that instrument, rounded to the nearest sample. We test the application of this additional delay both without (hereafter re-alignment, Figs A1d and e) and with frequency filtering of the geophone components (hereafter re-alignment + filter, Figs A1g and h) to equalize the signal content on the X and Y channels. Filtering is performed using a low pass zero-phase sine-squared filter with the upper roll-off from 150–300 Hz, which is above the bandwidth of interest for the analysis of the GI source SWS (Fig. A1f). For each of these three test cases: initial, re-alignment, and re-alignment + filter, tilt analysis and R-T rotation analysis was performed.

A.2 OBSERVATIONS

A.2.1 Instrument tilt

We observe that for OBSs where the delay value generated by the action of the X component functional filter is greater than 1 sample (0.25 ms), there is a decrease in the value of the standard deviation of the misfit of the direct arrival polarization azimuth to the shot-receiver azimuth. This arises since the calculation is highly sensitive to the windowing of the first half-cycle of the direct arrival, and, therefore, if this is misaligned, there is a systematic impact on the resulting calculation which results in a larger misfit. Across the data set, the range of overall tilt values is 1–8°, however of the 21 OBS which have well resolved values 16 are <5° and 9 are <3°. Two instruments, OBS 7 and OBS 9, do not produce valid solutions to the tilt analysis in the range ±10°, and are excluded from further analysis at this stage.

A.2.2 Rotation to R–T

For the three different test cases, initial, re-alignment and re-alignment + filter, we observe that for OBSs where the delay value generated by the action of the X geophone functional filter is greater than 1 sample (0.25 ms), there is an improvement in the R:T power ratio, indicating that more energy is being successfully removed from the T component, which may be by up to an order of magnitude. This improvement is seen for both the re-alignment and re-alignment + filter test cases. The effect of applying the corrections on the rotation angles is a tendency to reduce the amplitude (b coefficient) of the cos2θ trend observed with azimuthal variation. Any changes to the overall mean rotation angle (a coefficient) are very small, typically <1°. This observation suggests that overall the alignment of the peaks has a bigger impact than the differences in the frequency content.

For small tilt values, for example as seen in Figs A3(d)–(g) and A4(d)–(g), the differences between alike pairs of rotation tests, pre- and post-tilt correction, are small. This observation indicates that instrument tilts present in this dataset, of up to 8° but predominantly less than 5°, do not have a major impact on the fidelity of the OBS horizontal geophone records. However, we choose to apply the tilt corrections derived for each instrument to the data.

APPENDIX B: CALCULATION OF SHALLOW VELOCITY STRUCTURE AT SCANNER POCKMARK COMPLEX

B.1 INTRODUCTION

As discussed in Section 5, in order to relate observations of P–S converted arrivals (and shear wave splitting) in the GI data to observations of stratigraphic structure derived from P-wave reflection data, it is necessary to generate sufficiently finely depth-resolved models of subsurface P- and S-wave velocity. Here we describe the procedure used to derive P-wave velocity models for each instrument, which is then converted to S-wave velocity using the relationship of Castagna et al. (1985) for mudrocks.

B.2 METHOD

Instruments were selected for the analysis to give a range of geological settings, while also providing good resolution and simple solutions. Therefore, instruments and shot lines associated with subsurface dipping structures were excluded, in order to simplify the modelling stage to planar, horizontal, laterally homogeneous layers. Instruments located above shallow gas accumulations (Fig. 3) were also excluded, since the gas shadow effect means that reliable velocity values cannot be obtained for a time-interval below the gas. A combination of single, close-approaching shot lines and whole instrument gathers were used to generate semblance panels, with analysis offsets restricted to ±200 m. A total of nine instruments were used in velocity modelling, including OBSs 8, 12, 19 and 20 on which shear wave splitting is observed. OBS 1 was excluded due to the effects of the gas blanking reducing the ability to observe reflection events at intermediate depths. The instruments used in the analysis are shown in Fig. 2.

A sparker source is required for this analysis in order to achieve the sufficient resolution required to image the shallowest subsurface layers. However, Duraspark source streamer recordings show that this source cannot distinguish clearly between the UWG and LWG units (see Figs 3d–f), which is confirmed by the first major semblance peaks being located at ∼20 ms after the direct arrival on the OBS data. As a result, the Witch Ground is treated as a single layer in the velocity modelling. Normal moveout velocities were picked using a semblance approach, in order to flatten arrivals associated with P-wave reflection events recorded on the OBS hydrophone channel, from the Duraspark surface sparker source. A single composite NMO-time trend was picked for each instrument using a combination of close-approaching shot lines and the whole-instrument gather as appropriate (Figs B1a and B2a).

Figure B1.

Velocity analysis for Vp model construction, OBS 12. (a) Semblance panels for individual shot line (first and second panels) and whole instrument gather (third panel), offsets ±200 m. Black squares show location of composite velocity picks. (b) P-wave interval velocity–depth layered model constructed by conversion of moveout velocity using the Dix equation (1955). (c) Duraspark hydrophone record with normal moveout from (a) applied, showing multiple shot lines with closest approach <150 m, and total offsets shown ±500 m. Plotted in order of shot number within dataset. Red bars indicate locations of traveltimes picked on flat reflection events. (d) rayinvr modelled traveltimes (black) compared to picked traveltimes from data (red), plotted with normal moveout of picked data removed. Data in (a) and (c) are filtered with a 1–200–500–1000 Hz zero-phase bandpass filter.

Figure B2.

Velocity analysis for Vp model construction, OBS 19, plotted as for Fig. B1. In (a) panels 1–4 are individual shot lines and 5 is the whole instrument gather.

A layered modelling approach was chosen to test the resulting velocity models, since this would recreate the principles required to generate strong reflections in the data, by having velocity contrasts rather than simply changes in velocity gradient (Figs B1b and B2b). However, a reflectivity approach was not used for modelling. Within each model, the boundaries between layers are planar and horizontal and the layers themselves are laterally and vertically homogenous. For each OBS in the analysis, the subsurface velocity–depth model derived from the analysis above was used to generate the theoretical arrival times of reflections arriving from the base of each model layer using rayinvr (Zelt & Smith 1992). The model fit was determined by comparing the calculated traveltimes through the model with picked arrivals made from the data, using all available close-approaching shot lines, but with longer offsets, up to 500 m, permitted, in order to ensure enough arrivals can be picked (Figs B1c and B2c). Pick errors of 1 ms were applied uniformly throughout the data set, and the resulting χ² was calculated by comparing the times of the picked arrivals with through predicted by the ray tracing (Figs B1d and B2d). A χ² value of 1 corresponds to the misfit being equal to the assigned error of 1 ms. In general, we found that the models tend to over-fit the data, that is the χ² value is <1 (Figs B3 and B4).

Figure B3.

Model velocity sensitivity testing, OBS 12. Model layer numbers, corresponding to picked semblance peaks below seabed (model layer 1). Black squares show χ² fit values calculated using a pick error of 1 ms. Model layer numbers and number of picks are labelled above each panel. Note that the vertical scale for layers 1 and 2 is different. No perturbation is applied to layer 1, the water layer.

Figure B4.

Model velocity sensitivity testing, OBS 19, plotted as for Fig. B3.

To test each model's sensitivity to perturbation, we apply a scaling factor to the calculated model velocities, in the range 0.95–1.05, to give perturbations of up to ±5 per cent (Figs B3 and B4). There is a trade-off between the effects of varying model interface depths and velocities. For simplicity, and to reduce the total number of model runs, all layers within each model were perturbed equally in each test run. However, the χ² fits were calculated both for picks from each layer separately, and the model as a whole, in order to differentiate changes in sensitivity through the model.

B.3 OBSERVATIONS

We observe that the first semblance peak below the direct arrival, commonly arriving ∼20 ms after the direct arrival and corresponding to the base of the WG (see Fig. 3), often shows velocities with values less than the direct arrival moveout velocity, which is ∼1490 m s^–1, the velocity in water (Figs B1a and B2a). Below the WG the velocity rapidly increases, to values ≥1600 m s^–1, reaching a maximum of ∼1900–1950 m s^–1 in the AbG.

We observe that overall the model velocity sensitivity tends to decrease with depth, with perturbations of ±∼3–4 per cent required to cause significant divergence from the fit minimum for the WG, corresponding to model layer 2 (noting the different vertical scale of the panel). Perturbations of ±5 per cent result in χ² values of 1.66–2.27, corresponding to a rms misfit of 1.3–1.5 ms for OBS 12 (Fig. B3), while the fit is better for OBS 19 (Fig. B4). The perturbation required to produce significant misfit decreases in the LGM and CP (layer 3), although there are generally fewer picks which can be made in this region due to the less clearly internally layered structure of the sediments (Fig. 3). For layers corresponding to the LB and AbG, numbered greater than three, ±1–2 per cent perturbations are required to produce significant misfits, with χ² values ≥3, corresponding to an rms traveltime misfit of ∼2 ms, or double the error assigned to the picks (Figs B3 and B4).

In this analysis we concentrate on the effects of velocity uncertainty. To simplify the analysis, all layers within each model were perturbed equally in each test run. Thus, the effects of shallow layer perturbations do have an effect on deeper layers in the model. As a result, the perturbation limits determined from this process correspond to lower-bound estimates. We do not model the effects of depth uncertainty in our analysis, since the depth to individual layer locations is better characterized using the streamer reflections, rather than the time of the semblance peaks on the OBS data. We compensate for the effect in the P–S arrival time prediction step, where we average the individual models to produce a composite velocity–depth profile (Fig. 5), and then use this alongside picks from the sparker source streamer MCS records to get the best depth, and, hence, P–S time predictions of the horizons of interest. As a result of the averaging process and construction of a composite layered model, residual uncertainties in layer depths are reduced significantly.

Overall, we find that our analysis shows distinct variation of seismic velocity with depth, with the greatest increase being at the base of the WG. The shallowest layers are also the least well resolved, which is an impact of both the setting and the approach used. The moveout to interval velocity correction is highly sensitive to very thin layers, and thus small variations in the picked moveout can lead to relatively large variations in interval velocity. However due to the small path length, small velocity perturbations of ±3–4 per cent for the shallowest layer and ±1–2 per cent for deeper layers during sensitivity testing do not result in traveltime misfits which greatly exceed the picking error.