Geohazard detection using 3D seismic data to enhance offshore scientific drilling site selection

A geohazard assessment workflow is presented that maximizes the use of 3D seismic reflection data to improve the safety and success of offshore scientific drilling. This workflow has been implemented for International Ocean Discovery Program (IODP) Proposal 909 that aims to core seven sites with targets between 300 and 1000 m below seabed across the north-western Greenland continental shelf. This glaciated margin is a frontier petroleum province containing potential drilling hazards that must be avoided during drilling. Modern seismic interpretation techniques are used to identify, map and spatially analyse seismic features that may represent subsurface drilling hazards, such as seabed structures, faults, fluids and challenging lithologies. These hazards are compared against the spatial distribution of stratigraphic targets to guide site selection and minimize risk. The 3D seismic geohazard assessment specifically advanced the proposal by providing a more detailed and spatially extensive understanding of hazard distribution that was used to confidently select eight new site locations, abandon four others and fine-tune sites originally selected using 2D seismic data. Had several of the more challenging areas targeted by this proposal only been covered by 2D seismic data, it is likely that they would have been abandoned, restricting access to stratigraphic targets. The results informed the targeted location of an ultra-high-resolution 2D seismic survey by minimizing acquisition in unnecessary areas, saving valuable resources. With future IODP missions targeting similarly challenging frontier environments where 3D seismic data are available, this workflow provides a template for geohazard assessments that will enhance the success of future scientific drilling.

1 Introduction

When planning an offshore drilling campaign, one of the primary technical concerns that governs site selection is safety (Jeanjean et al., 2005; Mearns and Flin, 1995). The subsurface can be hazardous and is full of unknowns; therefore, it requires the full interrogation of all available data to reduce risks during drilling. Common subsurface geohazards include phenomena related to excess pore pressure such as shallow hydrocarbons, shallow water flows, faulting to shallow depths, mud volcanoes, and pockmarks, all of which can lead to incompetent sediments and seabed instability (Aird, 2010; Jensen and Cauquil, 2013; Ruppelt and West, 2004; Wood and Hamilton, 2002). Identifying these hazards prior to drilling allows for decisions to be made during site selection to either avoid the hazard completely, mitigate it or select the lowest risk option (Aird, 2010; Jensen and Cauquil, 2013). This process, often termed a geohazard assessment, is integral to drilling success and safety, with its effectiveness often relying on the availability of high-quality data coverage such as seismic reflection, bathymetry and well data.

The importance of geohazard identification is amplified in environments such as deep water, high pressure/high temperature, glaciated margins and frontier petroleum provinces where there is little prior drilling experience (Galavazi et al., 2006; Moore et al., 2007; Weimer and Pettingill, 2007). Within these locations, the frequency and significance of geohazards are often increased, as well as the consequences of accidents (Eriksen et al., 2014; Li et al., 2016). Furthermore, it is often these locations where data coverage and quality are poorest, resulting in a limited subsurface understanding and higher risk to drilling. In the last few decades, exploration for subsurface resources has expanded into more challenging environments due to increased global demand for petroleum products (Mitchell et al., 2012; Poppel, 2018; Suicmez, 2016). This has led to the acquisition of extensive 2D and 3D seismic reflection datasets and industry drilling in frontier areas. The expansion, especially with regards to new 3D seismic reflection data coverage, provides an opportunity to use these data for reasons beyond their original commercial purpose. This includes regional geological mapping and geohazard assessments allowing site identification and de-risking for scientific boreholes (Dutta et al., 2010; Selvage et al., 2012; Hovland et al., 1998).

A geohazard assessment involves the geospatial quantification of drilling hazards, ideally through the use of densely sampled 2D or 3D reflection seismic data (Aird, 2010; Khan et al., 2018; Selvage et al., 2012; Heggland et al., 1996). In the past, the required resolution to detect small but potentially hazardous features could only be provided through dedicated high-resolution site surveys across potential drill sites using such tools as sub-bottom profilers, small-volume airguns, or sparker reflection systems (Jensen and Cauquil, 2013; Parkinson, 2000). Today, the acquisition of large-scale 3D seismic reflection datasets in frontier basins provides spatial coverage that far exceeds most conventional site surveys (Games and Self, 2017). The processing of these datasets has improved sufficiently to provide a vertical resolution that approaches that of a traditional site survey and often provides a horizontal resolution exceeding that of even closely spaced 2D seismic site surveys (Games and Self, 2017; Oukili et al., 2019). Three-dimensional seismic surveys thus minimize the need for additional data, except for particularly complicated areas (Hill, 1996; Selvage et al., 2012; Sharp and Badalini, 2013; Williams and Andresen, 1996; Roberts et al., 1996).

Table 1 IODP expeditions and proposals with 3D seismic data.

Figure 1 Location map. A bathymetric map of the Melville Bay area in north-western Greenland from Newton et al. (2017), showing the extensive influence of glaciation on the seabed including the Melville Bay Trough (MBT) and Melville Bay Trough Mouth Fan (MB-TMF). The red square on the inset map shows the global location of the bathymetry. Annotations on the map include the distribution of regional rift elements, including sedimentary basins and structural ridges, as well as the location of paleo-shelf break positions of glacigenic prograding units from Knutz et al. (2019). IODP Proposal 909 site locations are shown as well as the associated stratigraphic targets of each site (I–VII). Shallow cores that were drilled by a consortium led by Shell that provide stratigraphic information used within the proposal are also shown (Acton, 2012; Nøhr-Hansen et al., 2018). The location of the key regional 2D seismic transect used for original site selection is also shown and represents the location of Fig. 2.

The increased availability of 3D seismic volumes in continental shelf areas often coincides with areas targeted by scientific drilling programmes. These include the International Ocean Discovery Program (IODP) that operates the drilling vessels Joides Resolution and Chikyu but also national facilities such as the Meeresboden-Bohrgerät (MeBO) of the Alfred Wegener Institute (AWI) and the British Geological Survey (BGS) Rock Drill. For this study, we present data from an Arctic frontier basin that support site selection for a proposed scientific drilling leg under the IODP. IODP has an excellent safety record, achieved through review of proposed sites by panels of international experts, and has drilled within frontier petroleum provinces in the past, with examples such as the early passive margin drilling along the western Atlantic margin (e.g. DSDP Legs 11 and 41; Ewing and Hollister, 1972; Lancelot and Seibold, 1977) and more recently in areas such as the Demerara Rise, offshore Suriname (Leg 207) and the Great Australian Bight (Leg 182) (National Research National Research Council, 2011). Several completed and proposed IODP expeditions report the availability of 3D seismic data within the study area, and it is possible that several others had undocumented access (Table 1). As commercial 3D seismic reflection datasets become more widely available for academic research, there is an opportunity to optimize both the scientific benefits and safety of proposed drilling sites by conducting more comprehensive geological and geohazard assessments. This study, whilst focussed on the hazards associated with IODP Proposal 909 within the north-western Greenland glaciated margin (Fig. 1), provides a geohazard assessment workflow that optimizes drill site selection through the use of 3D seismic data and serves as a template for the improved safety and success of future scientific drilling campaigns using 3D seismic data in frontier areas.

2 IODP Proposal 909

2.1 Setting

IODP Proposal 909 aims to drill a transect of seven sites across Melville Bay offshore north-western Greenland at water depths ranging from 0.5 to 1.9 km (Fig. 1). Here, a thick (> 2 km) Cenozoic sedimentary succession overlies a rift basin topography that formed during several stages of Early Cretaceous to Early Paleogene rifting between Greenland and Canada (Altenbernd et al., 2015; Gregersen et al., 2013, 2017; Oakey and Chalmers, 2012). This includes the extensive, elongate inversion structures of the Melville Bay and Kivioq ridges. These ridges separate the deep sedimentary basins of the Melville Bay Graben and Kivioq Basin and contain up to 9 km thick successions of syn-rift (seismic mega-units (mu) mu-G, -F and the lowermost -E) and post-rift (mu-E, -D, -C, -B and -A) sediments (Figs. 1 and 2) (Altenbernd et al., 2015; Gregersen et al., 2013, 2017; Knutz et al., 2015; Whittaker et al., 1997). Mu-E is attributed to the continental drift phase as seafloor spreading commenced in Baffin Bay. The lower part of mu-D is considered to have formed during the final syn-drift stage influenced by compressional tectonics as a consequence of Greenland converging with the North American Plate (Knutz et al., 2020). Deposition of the upper part of mu-D, representing a hemipelagic succession, was presumably deposited during a phase of post-drift tectonic relaxation. The mega-units are separated by seismic horizons (hz) e1, d1, c1 and b1 that are generally expressed as regional unconformities.

Figure 2 Regional geology. (a) A 2D regional seismic reflection line showing the stratigraphy and structure across the shelf. Regionally mapped seismic mega-unit interpretations from Gregersen et al. (2013, 2017) and Knutz et al. (2015) are shown as well as the projected locations of the seven IODP Proposal 909 stratigraphic targets and the comparative coverage of 3D seismic data. (b) An enlarged section of the seismic line from A that focusses on the glacigenic wedge that represents the Melville Bay Trough Mouth Fan system, showing the 11 interpreted sub-units of mega-unit (mu)-A from Knutz et al. (2019) as well as the projected location and depth of the IODP Proposal 909 sites. The location of the seismic line is shown in Fig. 1. (c) A stratigraphic column displaying seven seismic mega-units and their associated rift stages, as well as the expected age of the stratigraphy that comprises each of the stratigraphic targets (I–VII).

On the inner shelf margin, thick late Miocene and Pliocene (mu-C and mu-B) marine sediments constitute the uppermost post-rift sequence (Knutz et al., 2015). This includes widespread, late Neogene contourites and their correlative mass transport deposits, down-slope from a major erosional scarp above the Melville Bay Ridge (Figs. 1 and 2). The Neogene sediments are exposed at the seabed on the inner shelf and are progressively buried towards the basin by thick glacigenic packages forming part of mu-A (Fig. 2) (Knutz et al., 2019). The exposure and burial of the Neogene marine successions occurred due to multiple phases of ice sheet expansion since the late Pliocene. Through these glaciations, material was eroded and redistributed, leading to over ∼ 100 km of shelf edge progradation and accumulation of the Melville Bay Trough Mouth Fan (MB-TMF) (Figs. 1 and 2) (Knutz et al., 2019; Newton et al., 2017). These glacigenic progradational units likely consist of highly variable sediment lithologies and grain size compositions (Christ et al., 2020; Knutz et al., 2019). The prograding units are separated by unconformities that generate distinct unconformable seismic reflections that express detailed morphologies formed by sub-glacial erosion and deposition, similar to the features observed on the present seabed (Newton et al., 2020, 2017).

The glacigenic succession (mu-A) has been subdivided into progradational units that suggest a minimum of 11 major phases of ice advance and retreat across the shelf since ∼ 2.7 Ma. These sub-units (1–11), proposed by Knutz et al. (2019), are used throughout this study (Fig. 2b). Furthermore, the regional stratigraphic framework consisting of seven seismic mega-unit subdivisions (mu-G to -A) (Fig. 2) was proposed after extensive regional mapping of the north-western Greenland continental margin by Gregersen et al. (2013, 2017) and Knutz et al. (2015).

2.2 Scientific drilling objectives

IODP Proposal 909 aims to illuminate the late Cenozoic history of the northern Greenland Ice Sheet (GrIS) and specifically to ascertain the paleo-ice sheet dynamics during past warm climates (Knutz et al., 2018). This will be achieved by recovering drill cores at seven sites along a transect crossing the north-western Greenland margin (Figs. 1 and 2), each recovering key stratigraphic targets (Targets I–VII) to obtain a composite stratigraphic succession from Oligocene/Early Miocene to Holocene (Fig. 2). The overall objective is to examine the range of feedback and forcing mechanisms (oceanic, atmospheric, orbital, tectonic) impacting the GrIS through time – addressing several current themes of the IODP Science Plan (Bickle et al., 2011).

2.3 Site-selection requirements

The seven stratigraphic targets (Targets I–VII) were selected along the south-west–north-east trending regional 2D seismic transect (Figs. 1 and 2) and represent high accumulation rate deposits within the hemipelagic sequence of mu-D, contourite drifts of mu-B and -C, as well as potential interglacial and proximal shelf deposits within the trough mouth fan system of mu-A (Fig. 2). Once stratigraphic targets were defined, specific drill sites were selected. The strategy for site selection was to maximize both the chance of reaching the stratigraphic target and the chance of good core recovery whilst avoiding all identified potential drilling hazards. Additionally, the selection of several alternate sites was required in preparation for unexpected drilling issues or iceberg mitigation management.

2.4 Regional geohazard considerations

The north-western Greenland margin is a frontier petroleum province, with potential Cretaceous source rocks identified across the region in shallow cores and outcrop (Fig. 1) (Acton, 2012; Bojesen-Koefoed et al., 2004; Nøhr-Hansen et al., 2018; Núñez-Betelu, 1993) as well as deep unexplored sedimentary rift basins (Henriksen et al., 2009). The area experienced a surge in oil and gas exploration activity between 2007 and 2014, resulting in five exploration licenses being awarded within Melville Bay and the acquisition of extensive 2D and 3D seismic datasets. A shallow coring program was carried out in 2012, but no exploration wells have been drilled. Since then, the identification of a large potential gas reservoir (Cox et al., 2020a) and widespread evidence for shallow gas and gas hydrates (Cox et al., 2020b) further support an active petroleum system, thus underscoring the need for a geohazard assessment study prior to drilling.

Today, marine-terminating glacial outlets do not expand much beyond the coastline, and Melville Bay is generally free of sea ice cover during the summer, allowing access to industry and research vessels (Saini et al., 2020). However, Arctic weather conditions and the potential for icebergs carried northward by the West Greenland coastal current present a challenge to any logistical operations in the area. The environmental factors, considered alongside the possibility of shallow hydrocarbon occurrence and coarse-grained Quaternary sediments indurated by ice loading, makes drilling complicated and highly reliant on site survey data. To evaluate and select coring sites for Proposal 909, available 2D and 3D seismic reflection data were used to conduct a comprehensive geohazard assessment that maximized the chance of safe drilling and successful core recovery.

3 Data

The seismic data used within this study include four separate surveys that were acquired between 2007 and 2019 (Table 2). All data used were provided in SEG normal polarity with a downward increase in acoustic impedance represented by a red positive peak and a downward decrease in acoustic impedance represented by a blue negative trough (Fig. 2). The regional 2D data used in this study form part of four surveys that were acquired by the geophysical company TGS (2007–2010) as part of the Baffin Bay 2D regional dataset. The dataset was used to help understand the regional geology and map the spatial distribution of stratigraphic packages that represent drilling targets (full extent of the regional 2D survey shown by Gregersen et al., 2019 – Fig. 1). From the TGS 2D data, a subset of seven lines crossing the Melville Bay shelf was used directly within the site-selection process (Fig. 3, Table 2).

Table 2 Acquisition parameters and survey statistics for the four seismic surveys used within this study. A dash (–) represents a parameter that is not applicable for that survey. The Pitu HR parameters can be read from the Pitu survey. Abbreviations used within the table include two-dimensional (2D), three-dimensional (3D), ultra-high-resolution (UHR), reprocessed (repro.), number (No.), average (Avg) and two-way time (TWT). Units used include metres (m), kilometres (km), milliseconds (ms), seconds (s) and hertz (Hz).

^∗ The dominant wavelength was calculated using an average velocity of 2000 m s⁻¹.

Figure 3 A fluid hazard and data map of the study area displaying the regional bathymetry from Newton et al. (2017) as well as the location of all seismic data used within the Proposal 909 study. The locations of potential fluid anomalies were identified and mapped via the shallow gas detection process within the geohazard assessment workflow. Active and abandoned IODP Proposal 909 sites that exist within the 3D seismic extent are shown as well as the location of the Melville Bay Ridge (MBR) and Figs. 4a, b, 11a, b, 12e, 14a and b.

Two 3D seismic surveys represent the principle data used for geohazard detection within this study (Fig. 3). The first is the Pitu survey that was acquired by Cairn Energy PLC in 2011 (Table 2). This survey was also provided as a pre-stack depth migration (PSDM) volume. The survey was reprocessed by CGG in 2013 to provide the Pitu HR (high-resolution) survey, which is a subset of the full volume with increased spatial and vertical resolution (Table 2). The second 3D survey, the Anu survey, provides extensive 3D coverage towards the shelf edge and was acquired by Shell in 2013 (Fig. 3). In 2019, a new ultra-high-resolution (UHR) survey, the LAKO UHR survey (LAKO is an acronym for the research vessel HDMS Lauge Koch), was acquired across several preliminary site locations through collaboration between the Geological Survey of Denmark and Greenland and Aarhus University (Pearce et al., 2019) (Table 2 and Fig. 3). The aim was to provide UHR imaging of the upper 500 ms below the seabed in order to supplement the existing industry seismic data within the pre-defined target areas. In addition, several longer transects were obtained that could provide additional seismic–stratigraphic information and add to the reconnaissance of additional sites. Eight lines from the UHR survey were used to characterize drill sites located within the areas of 3D seismic coverage (Fig. 3).

Depth conversion

Depth conversion was required to provide accurate estimates of seismically defined drill target depths in metres. Velocity information was provided from industry wells drilled by Cairn Energy PLC ≳ 300 km south of the transect (location shown in Gregersen et al., 2019 – Fig. 1), from shallow core sites U0100/110 (Fig. 1) and from an interval velocity cube that was provided over the Pitu survey area (Figs. 3 and 4). The interval velocity cube was created from seismic velocities during seismic processing through the application of a Kirchhoff pre-stack TTI depth migration (PSDM) that used a bin size of 12.5 × 25 m and a migration half aperture of 4500 m (Fig. 4). This process produced a depth-converted version of the Pitu 3D survey. The interval velocities and the depth cube, however, have not been calibrated to measured depth data, due to large distances to the closest well (≳ 300 km), and likely contain some error.

Figure 4 Interval velocities created through Kirchhoff pre-stack TTI depth migration (PSDM) for the Pitu survey overlaid on seismic reflection dip (a) and strike lines (b) in two-way time across the Melville Bay Ridge structure. Velocity reduction likely due to the presence of gas can be observed in the location of the Eocene-aged reservoir on top of the ridge as well as beneath a bottom-simulating reflector (BSR) within the free gas zone. The locations of panels (a) and (b) are shown in Figs. 3 and 9.

Therefore, average velocities were determined using all of the velocity data available (mentioned above) and were then extrapolated across the study area using a comparison of the potential lithology and depositional setting of the sediments as well as their general depth. Linear time–depth equations were also generated for sites within the glacigenic wedge to consider elevated compaction and velocity due to ice loading. Velocities from this linear trend were then compared against the original estimated average velocities and interval velocities from the velocity cube to calculate the potential error (often ∼ 100 m). This error was used to adjust the final metric depth targets to minimize overestimation and avoid drilling into deeper, potentially hazardous intervals.

4 Geohazard assessment workflow

The following geohazard assessment workflow was used to select drill sites that represent the lowest possible risk whilst meeting the scientific objectives (Fig. 5). The workflow considers the increasing availability of data that are typical through a project's progression, starting with regional 2D seismic data to understand the regional geology and pick initial stratigraphic targets and sites. Three-dimensional seismic data are then used to conduct a more detailed interrogation of the subsurface in order to delimit geohazards within the proposal area and guide site selection and minimize risk. Finally, additional data, in this case 2D UHR seismic data, are used in collaboration with the 3D seismic to fine-tune the selected sites to ensure they represent the most suitable and safest locations possible.

Figure 5 A workflow diagram that outlines both the steps conducted within the geohazard assessment leading to site selection and how the differing seismic data types were used throughout the process. This workflow should be used as a guide for future geohazard assessments using 3D seismic data but is not fully exhaustive and can be amended to fit different datasets and locations through the addition of extra steps that fulfil the requirements of future projects.

4.1 Seabed

The first step of the 3D seismic geohazard assessment is to map the strong, positive amplitude event that represents the seabed across the area of 3D coverage (Fig. 6). This provided the seabed depth, firstly in TWT and subsequently in metric depth, after a conversion using a typical seawater velocity of 1480 m s⁻¹ (based on the Pitu interval velocity cube) (Fig. 6a). A bathymetric compilation from Newton et al. (2017) provided an additional high-resolution image of the seabed morphology. The mapping, supported by the bathymetry data, identified a wide range of seabed features that have been interpreted as being of glacial origin, created by recent shelf glaciations, including lineations, ridges, iceberg scours and near-circular depressions interpreted as iceberg pits (Fig. 6; Newton et al., 2017). These features can create localized areas of high-dip and highly compacted sediment and should be avoided to prevent instability of the coring equipment on the seabed (Bennett et al., 2014).

Figure 6 Seabed mapping. (a) A seabed structure map across the 3D seismic coverage shown in both two-way time and converted metric depth, with the bathymetry data from Newton et al. (2017) shown beneath. (b) A dip attribute map across the 3D seismic coverage that has been filtered to display features by dip severity, with colours amber and red relating to moderate (2–5^∘) and high (5+^∘) risk respectively. (c) An amplitude attribute extraction map at the seabed surface across the 3D seismic extent showing the variation in amplitude between areas dominated by or free from glacigenic features. The locations of both active and abandoned IODP Proposal 909 sites are shown in panels (a)–(c) along with the locations of panels (d) and (e), Figs. 7a and 13a. (d) A seismic cross section from the Anu survey showing an example of a glacigenic wedge which is topped by glacial lineations. (e) A seismic cross section from the Pitu HR survey showing the unconformable nature of the seabed as well as seabed depressions that likely represent iceberg pit marks (cf. Brown et al., 2017).

A structural dip attribute was extracted onto the mapped seabed surface to identify areas of high seabed dip (Fig. 6b). Filtering of the attribute allowed the severity of dip to be separated into areas of low (0–2^∘), medium (2–5^∘) and high (5^∘+) dip, which relates to areas of low, medium and high risk respectively. The seabed dip–risk cut-offs are used here primarily to help avoid glacial geomorphological features on the seabed, but this technique is applied in commercial geohazard assessments to consider the effect of the seabed structure on the critical failure of slope sediments as well as seafloor infrastructure tension and strength (Dan et al., 2014; Vanneste et al., 2014). The cut-offs selected are relatively low compared to other studies (e.g. Haneberg et al., 2015), but these were chosen to remain cautious and decrease potential risk, as much of the study area is relatively flat lying and therefore remains classified as low risk (Fig. 6b). Applying a cautious approach allows for steeper-dipping sections to be documented for further analysis should the underlying stratigraphy be potentially suitable for drilling – e.g. steep dip does not necessarily preclude drilling, but these bins allow for areas requiring greater consideration to be highlighted. Seismic amplitudes were also extracted onto the seabed surface (Fig. 6c). Extreme high or low amplitudes were shown to coincide with either glacial depositional features (such as Fig. 6d) (low amplitudes) or iceberg pit marks/depressions (high amplitudes) (Fig. 6e), supporting the inferred composition of elements represented by areas of high dip. Away from these features, the seabed amplitude is relatively uniform.

4.2 Faults

Fault mapping within the 3D survey area involved the use of the variance seismic volume attribute (coherency within other software packages) to image discontinuities, viewed mainly through time-slice intersections (z slice) (Fig. 7). Faults in the immediate area of the drill sites were manually mapped, including small faults belonging to a dense polygonal fault system within the post-rift stratigraphy of mu-D and -C (mainly within the Pitu area) (Cox et al., 2020a). Several deep-seated faults are observed extending close to the seabed (Fig. 7b), potentially connecting deeper fluid pressures to the shallow stratigraphy, which would represent a significant drilling hazard. Therefore, in areas close to a proposed site, fault characteristics were assessed based on their vertical extents, offsets and possible connections to deeper seismic anomalies (as shown in Fig. 7). High risks were associated with faults displaying connections to deeper anomalies; however, it was recommended to avoid all fault penetrations where possible to maximize the chance of good and stratigraphically continuous core recovery.

Figure 7 Fault mapping. (a) A time-slice intersection (z slice) through the Pitu HR survey displaying the variance structural attribute to highlight fault locations at a depth of −1258 ms two-way time. Interpreted fault planes (white dashed lines) are shown as well as the spatial distribution of a potential hydrocarbon fluid anomaly that may be in pressure communication to shallower depths due to the intersection with fault planes. IODP Proposal 909 sites MB-7A (abandoned) and MB-7B (primary) are also shown as well as the location of panel (b). The location of panel (a) is shown in Fig. 6a. (b) A seismic cross section from the Pitu HR survey in black and white to emphasize faults, which displays interpreted fault planes and the location of a deeper fluid anomaly as well as the location of panel (c). (c) An enlarged seismic cross section from the Pitu HR survey that shows the seismic character of the deeper fluid anomaly.

4.3 Shallow gas detection

Gas-related seismic anomalies are expected to display “bright”, anomalously high negative amplitudes due to gas fluids causing a reduction of the bulk modulus (due to the extremely low density of gas) causing a significantly negative acoustic impedance contrast across the boundary at the top of the gas-bearing reservoir (Cox et al., 2020c; Hilterman, 2001; Nanda, 2016). A positive seismic amplitude, often of comparable amplitude, is often associated with the gas–oil–water contact or the base of the reservoir if it contains anomalous fluids across its entire vertical extent. For thin reservoirs filled with gas, the response is often a highly asymmetric high negative to high positive doublet (Cox et al., 2020c; Raef et al., 2017). The presence of oil can cause a similar but often much reduced seismic anomaly. Due to these phenomena, the 3D seismic data were investigated to understand the seismic character of fluid-related anomalies within the study area. This included the physical character of anomalies, which were often isolated bright spots or brightening along single horizons (Fig. 8), as well as their spatial and stratigraphic distribution, which was widespread within all levels of the post-rift stratigraphy (mu-D to -A) and focussed around the Melville Bay Ridge (Cox et al., 2020b). The seismic amplitude range of the bright anomalies was also analysed to aid amplitude extractions. The analysis identified that an amplitude of negative (–) 10 000 or above likely represented a fluid anomaly, although some dim anomalies can occur beneath this limit, and this was considered during the filtering of results (Fig. 8c). Dimmer anomalies can be attributed to poorer reservoir quality or thin bed tuning effects, whereas amplitude is expected to be insensitive to gas saturation when this exceeds 10 % of the pore volume (e.g. Hilterman, 2001).

Figure 8 Shallow gas detection. (a) A 2D map view of several overlaid windowed minimum amplitude extractions from the Pitu HR survey that have been filtered to only display amplitudes that are greater than −10 000 and likely represent fluid anomalies. The extracted amplitudes have been colour coded to distinguish between windows from different mega-unit packages, giving a sense of anomaly depth, a technique which can also be applied to individual windows. The windows displayed are shown in both panels (a) and (b). The IODP Proposal 909 sites in this area are also shown in locations that avoid the potential fluid anomalies as well as the location of the seismic line used in panel (b). (b) A seismic cross section from the Pitu HR survey that highlights the main horizons used to create the proportionally distributed extraction windows. Several seismic anomalies that will be identified within the extractions are also shown along with the specific windows displayed in panel (a). (c) A histogram that shows the distribution of amplitudes from the selected extraction windows shown in panels (a) and (b). The histogram and the colour bars show how the amplitudes have been filtered to only show the extreme values which are more likely to represent fluid anomalies. These amplitude cut-offs can be altered for specific windows to either display or hide certain amplitude ranges.

Once this initial analysis was complete, the information was used to automatically extract amplitude anomalies from the 3D seismic data above the defined threshold that likely represents hydrocarbon occurrences. This process used a series of minimum amplitude extraction windows that cut proportionally through the stratigraphy, extracting the most negative amplitude at every seismic trace (Fig. 8). The amplitude extraction result was then filtered to remove the majority of data and highlight just the most negative amplitudes (those that likely represent fluid anomalies) (Fig. 8c).

The stratigraphy was split into two units for the extractions, the first between the seabed and seismic horizon (hz) d1 with 10 proportional extraction windows (each ∼ 50 ms TWT thick) and the second between hz d1 and hz e1 with 5 windows (each ∼ 35 ms TWT thick) (Fig. 8). Using relatively narrow windows provided some depth control for each anomaly, which can be displayed through colour coding (Fig. 8a), instead of having a wide depth estimate that is as thick as the seismic mega-unit. This was important when selecting a shallow drill site where many of the deeper anomalies well below the target are less relevant. It also allowed the windows to be displayed separately, i.e. just for a particular depth zone. Furthermore, the resulting amplitudes for each window were often filtered differently (often using opacity rendering) in an attempt to highlight the anomalies at that level and also to try and remove unwanted amplitudes that may mask fluid events, such as naturally bright horizons, data acquisition footprints (near the seabed) and amplitudes interpreted as not representing fluid anomalies.

Once the data had been filtered, the individual scrutiny of each anomaly was conducted to determine whether that anomaly was in fact fluid-related or whether the anomalous amplitude may instead represent something else (such as the high amplitude causing features mentioned above). This evaluation is subjective to some degree, and this should be considered in the final result. Nonetheless, this evaluation was important as it stops potentially credible locations from being ruled out due to non-fluid-related anomalies and may also highlight other dimmer amplitude anomalies that were missed. This whole process was iterative and involved returning to the attribute extraction or filtering stage to adjust the parameters when required to ensure all potential anomalies were considered.

Once the result was finalized, the locations of the remaining amplitudes were used both within the interpretation software directly and within GIS software for spatial analysis. The remaining amplitudes were either converted into 2D automatic boundary polygons or extracted as 3D geobodies to allow efficient visualization and to aid the creation of multi-layer hazard maps used within the site-selection process (Fig. 3).

4.4 Gas hydrate detection

Bright fluid anomalies were identified hosted within inclined, likely sandy, strata that terminate abruptly at a certain depth beneath the seabed (Fig. 9) (Cox et al., 2020b). This level, characterized by a dim, negative amplitude reflection that cuts across the stratigraphic layering in a manner that mirrors the seabed topography, represents a bottom-simulating reflector (BSR). Such cross-cutting features have been interpreted as representing the base of a zone containing stable gas hydrates, below which free gas can become trapped, resulting in a bright negative cross-cutting reflection that may be discontinuous or continuous depending on the degree of porosity and permeability variations of the host strata (Fig. 9) (Berndt et al., 2004; Kvenvolden, 1993).

Figure 9 Gas hydrate detection. (a) A composite seismic cross section through the Anu and Pitu HR surveys showing the distribution and seismic character of the identified BSR in the study area. The white dashed line (ca. 200 ms below seabed) at the BSR represents the pseudo BSR surface. (b) A 2D map view of a minimum amplitude extraction across the pseudo BSR surface showing the distribution of bright negative amplitudes that likely represent trapped free gas. The mapped area of the BSR is shown as well as the locations of the seismic lines used in panels (a), (c), and (d), Fig. 4a and b. (c) A seismic cross section from the Pitu HR survey confirming the presence of the BSR along with multiple free gas anomalies trapped beneath it. (d) A seismic cross section from the Pitu HR survey across the area containing a potential deeper BSR (ca. 300 ms below seabed). The location of panel (e) is also shown. (e) A zoom-in of a section of panel (d) showing the location of fluid anomalies (yellow lines) beneath the potential deeper BSR and much thicker gas hydrate stability zone (GHSZ).

BSR features were mapped within the 3D seismic surveys to identify an area of 537 km² that likely contains gas hydrates (Figs. 3 and 9). Bright, negative amplitude free gas anomalies are observed trapped beneath the BSR boundary and suggest a free gas column of up to 50 m thick (Fig. 9c). However, away from the bright free gas anomalies, mapping of the BSR became difficult due to an intermittent and dim BSR reflection. Therefore, due to the relationship between the BSR and the seabed, a pseudo BSR surface was created to aid the identification of this boundary in other areas of the survey (Fig. 9). This pseudo surface was based on empirical evidence from obvious BSRs in the data (Fig. 9c) and was created by cross-plotting the thickness of the GHSZ in areas containing these BSRs against the seabed depth. The best fit line between these two data was then used to create a surface at the expected BSR depth across the entire 3D survey extent, a surface that tolerates variations in stability zone thickness in response to changes in seabed depth.

This guide surface was used to focus BSR reconnaissance throughout the 3D seismic data. It was also used to identify additional areas of free gas (a much higher risk) that may be trapped at the base of the gas hydrates (Figs. 3 and 9). This involved combining the extensive pseudo surface with the BSR map and producing a minimum amplitude extraction for a 20 ms window across this surface (Fig. 9b). This technique, along with the manual reconnaissance of the seismic character at the pseudo surface depth, resulted in the identification of a previously undiscovered BSR as well as trapped free gas in the north of the Anu survey (Fig. 9a and b).

South of the Pitu survey, a cross-cutting negative amplitude reflection has also been interpreted as potentially representing a BSR; however, this feature exists ∼ 100 ms TWT deeper than the guide surface (Fig. 9a, d and e). This potential BSR is not as obvious as in other areas, but multiple stacked bright anticlinal anomalies exist within the underlying stratigraphy that may represent the upward flow of hydrocarbon fluids that are subsequently trapped beneath the BSR boundary (Fig. 9d and e). Rapid changes in BSR depth over short distances can occur and have been observed in other areas, e.g. the Lower Congo Basin (Andresen et al., 2011). However, as the seabed depth and shallow sediments are relatively consistent with the main BSR area a short distance away (< 5 km) (Fig. 9a), a significant variation in the phase boundary depth is unlikely. Therefore, these deeper anomalies could instead represent more traditionally trapped gas unrelated to hydrates (Fig. 9d and e).

4.5 Risking

Once the main geohazards described above were analysed, the individual results were combined to create a composite hazard map that was used directly for site selection (often known as a common risk segment – CRS – map) (Hill et al., 2015) (Fig. 10). The hazards considered within this process (seabed features, shallow gas and gas hydrates) were rated as either moderate risk, such as low dip glacial seabed features or areas containing gas hydrate (without free gas), or as high risk, such as seabed depressions or shallow gas occurrences. These features were then colour coded using a traffic light system, with red representing high risk, amber representing moderate risk, and green representing low risk. After initial site locations were proposed, significant faults were mapped in more detail around the site and subsequently added to the map to assess whether alterations were required.

Figure 10 A common risk segment (CRS) map showing the green, minimum risk, depth target zone for Target VII where sites can be located safely whilst meeting the stratigraphic objectives. This CRS map was used to guide the site selection and abandonment of the shown Target VII sites. Similar CRS maps were also created to guide site selection and abandonment for the other stratigraphic targets within 3D seismic coverage (Targets III–VI).

All areas surrounding the moderate and high risk features could be considered minimum risk (green), but an attempt was made to maintain a minimum radius from each site location of 500 m to the nearest identified hazard. For each stratigraphic target (I–VII), a target seismic horizon and the maximum tolerable drilling depth to reach that horizon were defined, creating a more localized target zone where the site had to be located. These target zones were overlaid on the CRS maps to create the green, minimum risk zone, which highlighted the area that could be drilled safely and still reach the desired stratigraphic target within the tolerated depth (Fig. 10).

4.6 Site selection: geohazards and lithology

The geohazard assessment and CRS map creation provided a minimum risk target zone to guide safe site selection. However, several other logistical and geological factors had to then be considered whilst fine-tuning the final location of each site.

This included an assessment and interpretation of the seismic character to give an idea of the potential lithology of the sediments (following work by several authors on lithology interpretations from seismic data, e.g. Sangree and Widmier, 1979; Frey-Martínez, 2010; Hilterman, 2001; Badley, 1985; Stewart and Stoker, 1990). An attempt was made to avoid certain seismic features, such as (1) chaotic packages of reflections, often down-slope from erosional scarps, that likely represent mass transport deposits; (2) structureless to chaotic, near-surface packages that may represent boulder-prone glacigenic tills or debris flows; (3) high, positive amplitude reflections at the top of glacigenic progradational units that may represent paleo-seabeds that have been indurated due to ice loading during the following glacial period; and (4) relatively high, negative amplitude reflections that are often wavy (within mu-B or -C) or laterally continuous within mu-D that may represent the top of sand packages. These seismically inferred lithologies were considered risks and avoided due to the likelihood of poor core recovery and the possibility of the coring equipment becoming stuck in the hole.

The stratigraphic dip (on seismic) was also considered in an attempt to try and target flat horizons to maximize the chance of good core recovery. Finally, the site priority (either primary or alternate) had an influence on location, as alternate sites attempting to reach the same stratigraphic target as the primary one were required to be sited a minimum of 5–10 km away from the next nearest site (primary or alternate) if possible to manage the iceberg risk. After considering these additional factors within the low-risk zone provided by the CRS map, a suitable drilling location was selected.

5 Assessment results – Proposal 909 sites

Seven primary and 15 alternate sites have been identified for the 909-Full2 proposal (Figs. 1 and 3) (Knutz et al., 2018). These sites cover the seven stratigraphic targets (Targets I–VII) that have been identified in order to meet the scientific objectives of the drilling proposal (Figs. 1 and 2). The selection of these sites was based on multiple data that have become available as the proposal developed since its inception in 2016. Initially, the regional 2D seismic data represented the principle source of information for identifying drilling targets along the key seismic transect (Fig. 2). Since then, the majority of site selections have been refined based on industry 3D seismic volumes and the LAKO UHR seismic data that were collected in 2019 (Table 2). In the final iteration of the proposal, none of the proposed sites exists in an area only covered by regional 2D seismic reflection data (Table 3). The final review resulted in several site alterations and additions as part of Proposal 909-Full2 to accommodate alternate site requirements for iceberg-prone waters. This version was subsequently accepted by the Scientific Evaluation Panel and by the Environmental Protection and Survey Panel of IODP.

Table 3 A site information table for the active and abandoned sites of IODP Proposal 909. Seismo-stratigraphic units are from Gregersen et al. (2013, 2017) and Knutz et al. (2015) and are shown in Fig. 2. The ticks indicate the seismic data coverage at each site, while the red coloured ticks denote which data were used primarily for initial site selection. The result of the geohazard assessment column signifies the impact the assessment had on that site and includes confirm (site approved after being initially selected on a form of 2D seismic data), amended (site was moved to this location from an abandoned site), select (site was chosen as an initial location) and abandon (site was deemed an unfit location). Abbreviations used include stratigraphy (Strat), sub-unit (su), primary (Prim), alternate (Alt), abandoned (Ab), seabed (sb), metres below sea surface (mbss), metres below seafloor (mbsf), ultra-high-resolution (UHR), geohazard assessment (GHA), lithology (lith) and total depth (TD). All site locations are shown in Fig. 1.

The coverage of the three seismic data types (2D, 3D and UHR) varies across the sites, and their availability restricts which data can be used within the site-selection process (e.g. Target I and II sites which are outside of 3D seismic coverage; Fig. 1 and Table 3). However, for all sites located within the 3D seismic extent, these data and the geohazard assessment were considered the principle data that either guided the selection of new sites or were used to amend and confirm the suitability of locations that were previously selected on 2D data (Table 3). This led to several sites being abandoned (Figs. 1 and 3 and Table 3). In cases where the UHR seismic data were used to select a new site (such as MB-17A or MB-7B), the spatial analysis of drill targets and hazards from the 3D seismic reflection data were still the primary method used for final site approval.

5.1 Targets I and II

The deep water sites for Targets I and II (ca. 1950–1800 m water depth) are located beyond the present-day shelf break and aim to recover a paleoceanographic record of a Pleistocene drift system associated with the MB-TMF (Fig. 1 and Table 3). Target I represents mu-A subunits 9, 10 and 11, with Target II comprising an expanded section of the stratigraphically underlying subunit 8 (Fig. 2). These two sites were initially selected using regional 2D seismic data, with the majority having been subsequently refined by LAKO UHR data, apart from site MB-23A, which was selected directly on the LAKO UHR data. For the remaining sites, an assessment for potential drilling hazards was conducted using the 2D datasets (regional and UHR) which confirmed their suitability, but as these sites exist outside of the 3D seismic coverage, the full geohazard assessment could not be applied during the site-selection process.

5.2 Targets III and IV

Target III sites on the southern flank of the Melville Bay Trough (ca. 500 m water depth) aim to recover potential glacial and interglacial intervals expected to be of Early–Middle Pleistocene age within top-set strata of the MB-TMF that onlap onto glacial unconformities within mu-A subunits 6, 7 and 8 (Figs. 1, 2, and 11 and Table 3). Target IV sites in the southern central part of the Melville Bay Trough (MBT) (ca. 600 m water depth) focus on similar top-set strata covering a stratigraphic interval of likely Early Pleistocene age corresponding to mu-A subunits 3–6.

Figure 11 Targets III and IV. Seismic cross sections from the LAKO UHR survey that highlight the glacigenic stratigraphy of Targets III (a) and IV (b) as well as the locations of potential fluid anomalies and sites MB-31A (primary), MB-8A, MB-30A (primary) and MB-4C (alternate). The locations of both panels (a) and (b) are shown in Fig. 3, and uninterpreted versions are provided in Fig. S3.

Three of the sites, including both primaries, have been selected predominantly using the LAKO UHR seismic, with the remaining two alternate sites using the regional 2D survey (Table 3). During this selection, potential hazards observed in the 2D surveys, including amplitude anomalies, were identified and avoided. All of the selected sites, however, exist within the 3D seismic coverage (Fig. 3) and therefore have been analysed as part of the geohazard assessment. Each site location and the depth of the target seismic horizon were assessed against each step of the geohazard assessment workflow through the process of CRS map creation. In this case, no site alterations were required, and all sites selected existed within the green, low-risk zone (Fig. 10), avoiding all identified hazards, for example the potential fluid anomalies shown in Fig. 11, by an acceptable radius (usually > 500 m).

5.3 Targets V and VI

Target V and VI sites, located in the central MBT (ca. 600–700 m water depth), aim to recover pre-glacial Neogene contourite drifts of presumed Early Pliocene age (Figs. 2, 3 and 12) and a limited portion of the overlying prograding sediments, which may reflect the earliest marine-based glaciations in north-western Greenland (mu-A subunit 1).

Figure 12 Target V. (a) A seismic cross section from the LAKO UHR survey showing the location of primary site MB-17A. Both depth options for this site are shown along with horizon interpretations (coloured) that were mapped throughout the area to allow the correlation of target stratigraphy with the other sites. The section of seismic that represents panel (a') is also shown, with panel (a') displayed beneath, next to panel (d), to show how both the seismic lines in panels (b) and (d) are connected to the seismic line of panel (a). (b) A LAKO UHR seismic cross section displaying the stratigraphic correlation between sites MB-17A and MB-6D. (c) A seismic cross section from the Anu survey shown for comparison between the two surveys. Fluid anomalies are shown and numbered and their appearance can be compared with the associated anomalies displayed in the LAKO UHR survey in both panels (a) and (b). The location of panel (c) is also shown in both panels (a) and (b). (d) A LAKO UHR seismic cross section displaying the stratigraphic correlation between sites MB-17A and MB-5B. (e) A location map for the seismic lines used in panels (a), (b) and (d) that represents an enlarged section of Fig. 3 displaying the regional bathymetry from Newton et al. (2017) as well as the location of seismic data and mapped seismic fluid anomalies (red). The locations of the seismic anomalies that are highlighted and numbered in panel (c) are also shown. The location of panel (e) is shown in Fig. 3. Uninterpreted versions of panels (a), (b) and (d) are provided in Fig. S2.

Analysis of both the new LAKO UHR data and the geohazard assessment from the 3D seismic reflection data led to the selection of a new primary site (MB-17A) that compared favourably to the originally proposed site (MB-5B). Site MB-17A fulfils the Target V criteria whilst avoiding several potential fluid anomalies and minimizes drilling through potentially boulder-prone glacial debris flow sediments (Fig. 12d). The new site also allows operational flexibility, providing an alternate site for Target VI by sampling both Target V and VI sediments (e.g. Option 1 and Option 2) (Fig. 12 and Table 3). The recovery of the youngest drift sediments has been optimized by targeting a sequence of evenly layered strata that are located less than ∼ 60 m below the seabed (Fig. 12a).

Figure 13 Target VI. (a) A minimum amplitude attribute extraction along a mapped horizon representing the top of the contourite mound package that constitutes Target VI (horizon shown in panels b and c). The attribute extraction highlights bright amplitudes that may represent either tuning or fluid anomalies. The spatial relationship between these bright anomalies and fluid anomalies extracted through the shallow gas detection process (red dashed polygons) (displayed in Fig. 3) is highlighted. The Target VI site locations are shown as well as the locations of panels (b) and (c). The location of panel (a) is shown in Fig. 6. (b) A seismic cross section from the Anu survey showing the contourite mound that comprises Target VI. The occurrences of the potential fluid anomalies are also displayed along with the location of abandoned site MB-6C. (c) A seismic cross section from the Anu survey showing the location of the Target VI primary site MB-6D. The site location is shown to penetrate the top mound horizon down dip from the potential fluid anomalies on the packages' flanks.

Target VI sites aim to recover the oldest stratigraphic section of the Neogene contourite drift within mu-B. For Target VI, in combination with Target V, the overall strategy is to obtain a composite, high-resolution record containing the Early Pliocene warm phase to the Late Pliocene cooling (Table 3). The main drilling target is an expanded section of the wavy-mounded contourite drift that accumulated over an underlying erosional unconformity (hz c1) (Figs. 2, 12 and 13). Both the primary and alternate sites were picked directly on the 3D seismic data using the geohazard assessment as a guide. Primary site MB-6D is also covered by LAKO UHR data which were used to confirm the site (Table 3).

Initially, the primary for Target VI was selected using regional 2D data but was located in an area just beyond the limit of the 3D seismic coverage (between the Anu and Pitu surveys – Fig. 3). Therefore, it was decided to relocate the site to within the 3D seismic extent to allow a more detailed site analysis using the geohazard assessment (to MB-6C – Fig. 13). The shallow gas detection analysis (Fig. 8) identified two bright spots at the edges of the mounded contourite target package that may represent tuning effects between the negative amplitude reflection at the target top and a short-extent, possibly cross-cutting, positive amplitude event beneath it, but it could also (in a less likely worst case scenario) represent pockets of gas-charged sand at the target top (Fig. 13b and c). MB-6C, however, targeted the central part of the mound that looked evenly stratified and free of bright events (Fig. 13b), but after further considerations of the amplitude distributions and target dip, the site was moved to a new position (MB-6D), where it would penetrate the potential gas-hosting sandy horizon further down dip. This was a cautious effort to reduce the chance of encountering gases that have migrated up dip (Fig. 13).

The geohazard assessment also highlighted a potential fluid anomaly (bright spot) directly in the area where MB-6D is located (∼ 180 m to the N/NE) but ∼ 50 m beneath the base of the targeted drift package (Figs. 3 and 13). Due to this feature, a conservatively shallow target depth was selected that maintains a depth stand-off of ∼ 80 m, minimizing the chance of drilling too deep and reaching the underlying anomaly (Figs. 12 and 13).

5.4 Target VII

Target VII sites are located in the inner central part of the MBT (ca. 750 m water depth) within the area covered by the Pitu 3D (and Pitu HR) seismic survey and comprise the deepest planned sites within the proposal (Fig. 1 and Table 3). Target VII represents an apparently continuous Miocene succession (including a middle Miocene unconformity – hz d1) that has been exhumed on the inner shelf (Figs. 3 and 14), down to the top of a sedimentary wedge of likely Oligocene age (the target horizon) (Figs. 10 and 14a) (Gregersen et al., 2019, 2013). All of the Target VII sites have been selected using the geohazard assessment (Fig. 10 and Table 3).

Figure 14 Target VII. (a) A seismic cross section from the Pitu HR survey that displays the location of primary site MB-7B, the alternate sites MB-16A, MB-12A and MB-11A and the abandoned site MB-7A. Regional seismic mega-unit interpretations are shown as well as the target seismic horizon used to define the depth target used within the CRS map in Fig. 10 that represents the top of a potentially sandy MTD package. (b) A seismic cross section from the LAKO UHR seismic that shows potentially hazardous seismic features near site MB-7A that are also observed in the Pitu HR survey (shown in panel a) and ultimately led to the abandonment of that site. A projection of primary site MB-7B is also shown penetrating a section containing continuous and horizontal seismic reflections. The locations of panels (a) and (b) are shown in Fig. 3, and uninterpreted versions are provided in Fig. S3.

The previous primary site, MB-7A, had been selected on the key transect of regional 2D seismic (Figs. 2 and 14). However, the geohazard assessment, in combination with the LAKO UHR seismic, identified a number of issues with the original site location that were not identified using the regional 2D seismic.

1. Narrow, vertical sections of acoustic blanking were identified in both the 3D and LAKO UHR seismics that have been interpreted as possibly representing fluid-flow pipes (cf. Cartwright et al., 2007) or gas streaking (Fig. 14). Although site MB-7A did not intersect these features, they were within close proximity (< 500 m).

2. The MB-7A location was within close proximity to a deep-seated fault that potentially could create pressure communication between the shallow subsurface and a deeper, anticlinal anomaly that may represent trapped hydrocarbon fluids (Fig. 7).

3. The seismic signal is locally disturbed around the MB-7A site in both the 3D and LAKO UHR seismics (Figs. 5a and 7). This is likely caused by an overlying chaotic package of sediments directly beneath the seabed, affecting the signal beneath and possibly containing heterogeneities that would affect drilling, but this may also suggest a higher sand content at the target interval.

4. Potential thick mass transport deposit (MTD) sand packages were identified in the lowermost section of the MB-7A site, both beneath the highlighted target horizon and onlapping onto it (Fig. 14). These sediments would likely cause poor core recovery and could contain high fluid pressures due to their structural dip and proximity to the ridge flanks (Figs. 1 and 14). The top of the MTD sand package was used as the target seismic horizon within the geohazard assessment and the site-selection strategy was to avoid drilling into these sands, maintaining realistic maximum drill depths and including as much of the overlying Miocene section as possible (Figs. 10 and 14).

These issues led to site MB-7A being abandoned and relocated to the new primary site MB-7B, which targets an area of flat-lying strata within a 1–1.5 km wide fault block (with no fault intersections) that is bounded by a combination of north-west–south-east trending deep-seated faults and the polygonal fault system (Figs. 7, 10 and 14).

6 Discussion

A geohazard assessment is conducted prior to drilling in an attempt to restrict unnecessary delays, reduce costs, avoid poor data collection and, most importantly, reduce the likelihood of dangerous drilling events (Aird, 2010; Nadim and Kvalstad, 2007). The assessment allows the selection of sites that represent the lowest possible risk whilst also achieving the scientific objectives. This requires a detailed spatial analysis of all potential risks and the consideration of additional viable target areas, both regionally and stratigraphically (Selvage et al., 2012). For IODP Proposal 909, the sensitive environment associated with high-latitude continental shelves, as well as the likelihood of hydrocarbon occurrences, made a robust risk analysis increasingly important (Hasle et al., 2009; Nadim and Kvalstad, 2007; Li et al., 2016). The geohazard assessment was conducted in line with commercial site safety analyses (Jensen and Cauquil, 2013) whilst focussing on hazards that commonly create risks to drilling operations within both deep water continental margin settings (hydrocarbon occurrences, gas hydrates, near-surface faults, etc.; Aird, 2010; Jensen and Cauquil, 2013; Minshull et al., 2020; Ruppelt and West, 2004) and glaciated margins (glacial seabed features, problematic lithologies, indurated horizons, etc.; Bennett et al., 2014; Newton et al., 2017).

The assessment identified a pervasive distribution of hydrocarbon-related anomalies across the 3D seismic coverage which exist within all levels of the post-rift stratigraphy (mu-D to -A) (Figs. 3 and 8–10). The majority of these anomalies exist within the shallow subsurface (top 1 km of sediment) and most likely represent pockets of trapped gas or gas hydrates (Hilterman, 2001; Nanda, 2016). It was imperative to identify and avoid the shallow fluid anomalies prior to drilling as unexpected pressure kicks caused by low-density hydrocarbons can often lead to shut-ins and site abandonment and, in a worst case scenario, can cause blow-outs to occur (Holland, 1997; Prince, 1990). Gas hydrates, however, can be drilled through successfully and have even been the focus of several recent coring campaigns (Khabibullin et al., 2011; Ruppel et al., 2008; Ruppel, 2018; Wei et al., 2019). Though the hydrate deposits in this area are underlain by a free gas column that is up to 50 m thick (Fig. 9) (Cox et al., 2020b) and although overpressure beneath the hydrates is unlikely (thought to be at hydrostatic pressure; Dickens and Quinby-Hunt, 1994), the free gas content may create significant buoyancy pressure which could be catastrophic if drilled (Holland, 1997; Minshull et al., 2020). Therefore, all potential hydrocarbon fluid anomalies were classified as “maximum risk” and avoided as a priority (Fig. 10).

A dense network of both tectonic and polygonal faults was also identified within the study area (Figs. 2, 7 and 10). Avoiding fault intersections while drilling is important as deep-seated near-surface tectonic faults can pass high fluid pressures and hydrocarbons along the fault plane and, if not controlled, can lead to a loss of borehole control (Fig. 7) (Frydman et al., 2017; Jensen and Cauquil, 2013). Fault zones also have a much lower fracture gradient than non-faulted zones, which may lead to sediments breaking up during drilling and poor core recovery. Polygonal faults, however, are thought to be sealing under non-extreme pressure conditions (caused by factors such as glacial loading) and exist within many effective petroleum system seals worldwide, thereby negating the risk of fluid flow along the fault plane (Cartwright, 2019). Still, densely polygonally faulted successions close to the seabed can hold significant fluid pressures, although it is likely that the shallow fluid system here would have been depressurized during repeated glacial loading and unloading cycles (Fjeldskaar and Amantov, 2018; Goulty, 2008; Ostanin et al., 2017) that occurred from 2.7 Ma (Knutz et al., 2019). Polygonally faulted clays with high smectite content can be problematic, however, as they can cause drill pipes to become stuck due to the clays swelling when in contact with borehole fluids (Anderson et al., 2010).

Lithology variations were identified during the fine-tuning of site locations in an attempt to avoid problematic lithologies such as glacial debris flow deposits which may contain boulders (potentially damaging coring equipment), hard indurated horizons or mass transport deposits which affect the continuous chronostratigraphic nature of successions due to sediment re-deposition (Bennett et al., 2014; Jensen and Cauquil, 2013). An attempt was also made to avoid sandy horizons that were identified through certain acoustic impedance contrasts (often semi-bright, negative amplitude events denote a top sand horizon) or through a chaotic seismic facies (e.g. Fig. 14). These uncompacted horizons would likely cause poor core recovery, possibly exacerbated by high fluid pressures within the sands which would further promote sediment collapse (shallow water flow) (Ruppelt and West, 2004). Uncompacted coarse sand could even lead to the coring equipment getting stuck, such as was experienced on ODP Leg 174A (Austin et al., 1998).

Once a potential hazard was identified, the primary concern was to eliminate any drilling through or near to the potential hazard, more so than to unravel the detailed nature of the feature observed. However, the detailed results provided by the geohazard assessment coupled with a desire to understand the complex geological and fluid migration history of the area allowed a closer assessment of the identified hazards, with the most notable conclusions being discussed below.

Figure 15 A schematic summary of the structure, stratigraphy, geohazards, and fluid flow history that has been analysed and used within the process of site selection for IODP Proposal 909. Key features related to the elements of the geohazard assessment workflow are shown, including the location of hydrocarbon fluids, sandy lithologies, faults, and the potential source rock horizons and possible fluid migration pathways. The surrounding thumbnails show how these features were analysed and considered within the geohazard assessment workflow. Highlighting the distribution of these features across the schematic aims to provide context on how these features fit together and interact within the subsurface. The understanding of this system, along with the distribution of potential geohazards, was crucial in the identification of viable stratigraphic targets and safe drill sites.

Firstly, the varying seismic character of the identified fluid anomalies (Figs. 3 and 8) along with their relationship with stratigraphic and structural elements (Fig. 15) suggests various styles of trapping mechanisms along with a complex fluid migration history (Figs. 2, 8, 9 and 15) (Cox et al., 2020b). The distribution of fluid anomalies shows a concentration above the Melville Bay Ridge, suggesting that this positive relief feature focussed the upward migration of hydrocarbons in this area (Figs. 2, 3 and 15). This conclusion is further evidenced by the discovery of an extensive, likely gas-charged, Eocene-aged reservoir on the crest of the ridge (Figs. 8 and 15) (Cox et al., 2020a). Cretaceous source rocks are expected to exist within the buried syn-rift stratigraphy of the Melville Bay Graben and Kivioq Basin (Bojesen-Koefoed, 2011; Gregersen et al., 2019; Planke et al., 2009), which onlap onto the ridge high, thus providing the potential for hydrocarbon migration along sandy carrier beds or fault planes either up towards the ridge crest (charging the Eocene reservoir) or into the shallower post-rift stratigraphy (Figs. 2 and 15). Furthermore, the stratigraphy directly above the Eocene reservoir displays evidence for gas leakage, which is likely linked to the onset of multiple cycles of glacial loading and unloading of the crust and the mass redistribution of sediment associated with the development of the MB-TMF (Cox et al., 2020a, b; Knutz et al., 2019; Newton et al., 2017). These processes likely caused episodic variations in subsurface conditions and structural tilt that may have promoted pulses of fluid leakage (Figs. 1, 2 and 15) (Fjeldskaar and Amantov, 2018).

These processes also contributed to the concentration of hydrocarbons within the stratigraphy directly overlying the Melville Bay Ridge, causing the majority of this zone to be unsuitable for scientific drilling (Fig. 3). Moreover, expected fluid migration pathways from the deep basins and leakage from the Eocene reservoir are likely connected to the presence of gas hydrates, again within areas directly overlying the Melville Bay Ridge (Figs. 2, 3, 9 and 15). It is likely that the hydrate-forming fluids followed similar migration pathways to that of the identified gas anomalies, with free gas anomalies observed trapped at the base of the GHSZ at the present day, suggesting continued hydrocarbon migration and post-hydrate formation (Figs. 9 and 15) (Cox et al., 2020b). A portion of these fluids is possibly sourced from the leaking Eocene reservoir, but the presence of free gas anomalies within regional sandy horizons suggests that upward migration of these fluids may have occurred away from the ridge before entering these horizons and migrating laterally up dip into areas uplifted above the ridge structure (Figs. 9 and 15). Numerous fluid anomalies also exist in areas away from the ridge (Fig. 3), suggesting a more complex regional migration history which is likely characterized by the trapping and subsequent remigration of fluids in areas above the deeper basins (Cox et al., 2020b; Grecula et al., 2018). In addition, biogenic gas production, possibly from organic horizons within the contourite succession (mu-B and -C), may represent the source of gas observed within glacigenic progradational sands (mu-A) (Figs. 3 and 9b) (Muller et al., 2018; Rebesco et al., 2014). Anomalies away from the ridge, however, are more sporadic, and there are no identified gas hydrates or deeper gas reservoirs within the section investigated (Fig. 15).

Although not part of the principle task of identifying geohazards to minimize risk, the conclusions drawn from the more detailed assessment of the nature of the geohazard features identified provided a greater understanding of fluid distribution and an enhanced prediction of what identified seismic features may represent. This knowledge ultimately informed the assessment of risk and actively affected site-selection decisions and stratigraphic target amendments by avoiding areas characterized as representing focussed pathways for historical fluid migration (Fig. 15).

6.1 Influence on IODP Proposal 909

On the north-western Greenland continental margin, the geohazard assessment workflow was used to identify and document a wide range of potential drilling risks (Figs. 3, 5–10, and 15), which were considered on par with the scientific objectives and expedition logistics (number of drilling days available) when defining the final site localities. This process led to an efficient and informed selection of the primary and alternate sites, and the added detail and understanding provided through the assessment positively influenced the progression of the proposal through several stages of review and was integral to its success when considering its location within such a challenging region.

The influence of the assessment on minimizing risk is highlighted through a re-assessment of several original sites that were selected using regional 2D seismic data against the spatial distribution of geohazards identified through the 3D seismic assessment (Figs. 3, 7, 10, and 13–15 and Table 3). An attempt was already made to avoid potential drilling risks when selecting sites on the 2D data, but the restricted spatial coverage and lower seismic resolution led to a number of features remaining unseen in the 2D assessment, such as close-proximity out-of-plane fluid anomalies and faults (Figs. 3 and 7) or variations in seismic character suggesting sandy lithologies (Figs. 12 and 14). These types of geohazards were identified through the 3D assessment (and augmented by the LAKO UHR survey), leading to the alteration of these original sites to new and likely safer locations (Figs. 12–14 and Table 3). These alterations highlight the benefit of using 3D seismic analysis to identify safer areas outside of the 2D seismic lines whilst achieving the same or more optimal target parameters. Most importantly, the geohazard analysis enabled refined site selection within the Target VII area, which would have likely been abandoned if only covered by 2D seismic due to the deep drilling target, shallow gas distribution and a dense polygonal fault network, coupled with deeper faults and deeper gas sands (Figs. 2, 7, 10, 14 and 15).

The results of the 3D geohazard assessment workflow (along with the location of the regional 2D seismic) also significantly influenced the acquisition plan for the 2D UHR seismic survey. This survey was acquired in 2019 to complement the 3D assessment by providing a higher-resolution image of the top ∼ 500 ms below the seabed (Figs. 3, 11, 12, 14 and Table 2). These additional data both confirmed features observed on 3D seismic data and highlighted additional, more subtle features (such as potential lithology changes) which prompted several site amendments. Importantly, tailoring the acquisition using the 3D assessment results allowed the UHR survey to be focussed over areas that contained either the selected sites (to confirm location), more complex geology or a concentration of geohazards and also for target areas outside of the existing 3D seismic coverage (Targets I and II – Fig. 1). This allowed a more efficient integration of the two data types whilst minimizing acquisition within unnecessary areas, saving valuable expedition time.

6.2 Benefits to future projects

The geohazard assessment workflow provides a blueprint that maintains the high level of safety assessment associated with IODP drilling by incorporating modern 3D seismic data manipulation and interpretation techniques to take full advantage of the available data (e.g. Heggland et al., 1996; Sharp and Samuel, 2004; Sharp and Badalini, 2013) (Fig. 5 and Table 1). This helps improve both safety and the chance of success of future scientific drilling. Improved 3D seismic acquisition and processing workflows have led to significant improvements in 3D resolution and a reduced reliance on dedicated site surveys (Games and Self, 2017; Oukili et al., 2019). Comparison of the 3D seismic volumes as pseudo site surveys within Proposal 909, against the recently acquired 2D UHR site survey lines, shows that the quality and resolution of imaging are mostly consistent across the two data types, confirming its suitability (comparison in Figs. 12 and 14 and Table 2). The main differences observed included subtle seismic variations, possibly due to lithology differences noted on the UHR seismic, as well as some additional brightening observed in the Anu 3D survey, possibly caused by tuning effects from thin beds that are thicker than the resolution threshold on the UHR seismic and therefore do not display as brightly (seismic survey vertical resolution shown in Table 2 and Fig. 12b and c; Francis, 2015; Marzec and Pietsch, 2012). The 3D seismic volumes also meet the minimum acceptability criteria set out by the International Association of Oil & Gas Producers (IOGP) for using exploration 3D seismic data for site survey studies (Jensen and Cauquil, 2013).

Still, many projects never gain access to 3D seismic data, but the techniques defined within the workflow are still applicable in the 2D sense, especially when using high-resolution and closely spaced 2D site survey data instead. Therefore, the seismic interpretation techniques presented could help maximize the use of available data whilst increasing interpretation quality and efficiency for a wide range of future projects (Selvage et al., 2012; Smith et al., 2019). These methods include using amplitude extraction windows to quickly identify potential fluid anomalies (Fig. 8), using pseudo BSR surfaces created using the BSR–seabed relationship to expand and guide BSR interpretations (Fig. 9) or using stratigraphic depth target zones within CRS maps to collate all identified hazards while focussing site selection within an appropriate, low-risk area (Figs. 10 and 14).

The workflow also outlines how the increasing availability of seismic data through a project's progression can be combined to further assess potential geohazards and provide increased confidence in the chosen sites (Fig. 5). Importantly, the workflow addresses the key hazards that are likely considered within all site-selection safety assessments (Figs. 5 and 15), though the list is not fully exhaustive and several other geohazards or geophysical techniques could be included dependent on context and whether additional data are available. Additional datasets could include acoustic or elastic impedance volumes and derivative rock property estimations (e.g. Dutta, 2002; Huuse and Feary, 2005), “fluid” volumes, electromagnetic-derived resistivity volumes (e.g. Weitemeyer et al., 2006), angle stacks for use in amplitude versus offset/angle (AVO/AVA) studies (Castagna and Swan, 1997), etc.

Examples of additional or modified techniques include using stratigraphical (instead of proportional) attribute extraction windows to follow reflection dip and remove interference from bright amplitude regional reflections and also predicting the BSR depth from the water depth in areas conducive to gas hydrate formation but where no obvious BSR is observed (Field and Kvenvolden, 1985; Kvenvolden et al., 1993; Gehrmann et al., 2009). Alterations to selected cut-offs (such as seabed dip and amplitude filters; Figs. 6 and 8) and other parameters used to define the site location (such as the minimum radius to the nearest hazard) within the geohazard assessment workflow can also be made. If required, this may even include expanding the area where a specific site can be drilled, as seen used in regions even more prone to sea ice and icebergs, such as for IODP Exp. 379 (Amundsen Sea) (Gohl et al., 2019), where the safety radius (in this case along the 2D seismic line) defined around each site locality as being free of hazards was used as a zone of which last minute borehole adjustments could be made in response to the predicted ice trajectory.

The Gohl et al. (2019) study provides an example in which the subsurface- or geohazard-related requirements of a particular location would have influenced the techniques and processes considered within the workflow presented here. Other location-specific requirements include areas containing deep and active fault zones (e.g. the Nankai Trough, Japan) or hydrocarbon fluid venting and extensive gas hydrate deposits (e.g. the Cascadia Margin, offshore Oregon, USA). Drilling safely within such locations with specific geohazard requirements may only be possible when 3D seismic data are available for the assessment of geohazards and drill site selection. This is why IODP legs targeting such areas (Legs 322 and 311; Table 1) have been supported by dedicated, research council funded, 3D seismic acquisition and processing (Bangs et al., 2009; Scherwath et al., 2006). This again demonstrates the importance of 3D seismic data analysis for the future of scientific drilling.

7 Conclusions

IODP Proposal 909 aims to drill a transect of seven sites across the north-western Greenland continental margin that represents an area that is both a frontier petroleum province and a glaciated margin. A geohazard assessment was conducted that was optimized by the use of high-resolution 3D seismic data to accurately extract, document, and spatially analyse potential geohazard evidence and select drill sites that represent the lowest risk possible whilst meeting the scientific objectives in terms of realistic drill depths, section age, completeness, and thickness. The workflow undertaken for this assessment used 3D seismic analytical techniques to identify geohazards such as seabed features, fluid anomalies, faults, and certain lithologies. The mapped geohazards were combined to create common risk segment maps for each of the primary site groups, using a restricted (green) zone defined by the depth to the stratigraphic target to focus site selection. This process led to the alteration or abandonment of several sites that were originally sited on regional 2D seismic data to more optimal locations within the 3D seismic coverage. The workflow and results were also used to tailor the acquisition of an ultra-high-resolution site survey which optimized acquisition and overall expedition planning. This survey was then used in combination with 3D seismic to verify site locations and identify more subtle shallow features such as potential lithology changes. Ultimately, the full geohazard assessment workflow was used to support 7 primary and 15 alternate sites for the scientific drilling proposal.

The detailed, accurate and comprehensive results provided by the 3D geohazard assessment, as well as its influence on the success of IODP Proposal 909, highlight the importance and benefit of maximizing the use of all available data during the planning of a drilling campaign within a challenging environment. With both past and future IODP campaigns targeting areas such as frontier petroleum provinces or glaciated margins whilst having access to 3D seismic data, the detailed consideration of geohazards outlined by this workflow provides a template for future projects. Incorporation of 3D seismic data analysis into site selection and hazard evaluation will help allow a more comprehensive safety assessment that could enable scientific drilling in areas otherwise considered too risky whilst maintaining the high safety standards required by the IODP.