Shallow gas and gas hydrate occurrences on the northwest Greenland shelf margin

An extensive 3D seismic dataset was used to investigate the contemporary hydrocarbon distribution and historical fluid migration in Melville Bay offshore northwest Greenland, providing the first inventory of shallow gas and gas hydrate along this part of the Greenland margin. The shallow gas anomalies vary in seismic character and have been subdivided into four categories that represent (I) isolated shallow gas, (II) free gas trapped at the base of the gas hydrate stability zone (GHSZ), (III) gas charged glacial clinoforms and (IV) a giant mass transport deposit gas reservoir. Gas hydrate deposits have been identified across an area of 537 km2 via the identification of a discontinuous bottom simulating reflector (BSR) that marks the base of the GHSZ. The BSR has been used to estimate a geothermal gradient of 49 C/km across the GHSZ and a heat flow of 70–90 mW/m2, providing the first publically available heat flow estimates offshore western Greenland. The contemporary hydrocarbon distribution and historical fluid migration is influenced by the underlying paleo-rift topography and multiple shelf edge glaciations since ~2.7 Ma. Continued uplift of the Melville Bay Ridge, as well as glacial-sediment redistribution and basinward margin tilting from isostatic compensation, have led to a concentration of gas within the Cenozoic stratigraphy above the ridge. Furthermore, repeated variations in subsurface conditions during glacialinterglacial cycles likely promoted fluid remigration, and possibly contributed to reservoir leakage and increased fluid migration through faults. The top of the gas hydrate occurrence at 650 m water depth is well below the hydrate-free gas phase boundary (~350 m) for the present bottom-water temperature of 1.5 C, suggesting this hydrate province mainly adjusted to glacial-interglacial changes by expansion and dissociation at its base and is relatively inert to current levels of global warming. Glacial-related dissociation may have significantly contributed to the numerous free gas accumulations observed below the GHSZ at present day.

For these reasons, it is important to document and understand the distribution of shallow gas and gas hydrates, especially within high latitude environments sensitive to environmental change (Portnov et al., 2016). Understanding high latitude hydrocarbon systems can be difficult, however, as glaciations during the Pliocene-Pleistocene may have significantly influenced fluid migration (Goffey et al., 2016;Medvedev et al., 2019). The redistribution of sediment across glaciated margins, along with repeated ice loading and unloading on the shelf through glacial-interglacial cycles, can cause substantial structural changes, isostatic compensation, and significant variations in subsurface conditions (Fjeldskaar and Amantov, 2018;Zieba and Grover, 2016). These changes often promote the remigration of hydrocarbons, through processes such as fault reactivation, trap spill and seal breach, potentially leading to expulsion at the seabed into the water column (Goffey et al., 2016;Kjemperud and Fjeldskaar, 1992;Ostanin et al., 2017). Gas hydrate stability is also affected by variations in subsurface temperature and pressure (Kvenvolden, 1993), such as those related to glacial (un) loading, and can possibly lead to dissociation of hydrate at the seabed Grassmann et al., 2010). Attempting to understand this cryosphere-methane interaction may reveal the sensitivity of gas hydrate deposits to environmental change; providing critical insight into how these deposits may respond to future oceanic warming (Biastoch et al., 2011;Krey et al., 2009;Ruppel and Kessler, 2017).
The Melville Bay continental margin, offshore northwest Greenland, has experienced multiple episodes of ice sheet advance and retreat across the shelf since ~2.7 Ma, with sediment redistribution contributing to ~100 km of shelf edge progradation (Figs. 1 and 2) Newton et al., 2019Newton et al., , 2017. The region is thought to contain hydrocarbon systems (Bojesen-Koefoed, 2011;Gregersen et al., 2013;Whittaker et al., 1997) that were likely impacted by these glaciations. Here, extensive 3D seismic reflection data acquired for petroleum exploration are used to document the structure, stratigraphy and hydrocarbon occurrence across the glaciated margin of Melville Bay. The data coverage and mapping provides the most widespread and highresolution imaging of the subsurface and hydrocarbon anomalies along the west coast of Greenland, as well as the most northerly documentation of gas hydrates in Baffin Bay (Minshull et al., 2020;Nielsen et al., 2014). The primary objectives were to address: (1) Have hydrocarbons accumulated within shallow, post-rift sediments on the Melville Bay shelf; (2) What does the seismic character and distribution of observed hydrocarbon occurrences tell us about migration history; and (3) How have these accumulations been affected by recent glaciations? In addition to answering these questions, the analysis contributes to the wider understanding of Arctic gas hydrates, the West Greenland subsurface thermal regime and shallow subsurface drilling hazards for the upcoming IODP Leg (909) .

Regional Setting
Baffin Bay represents a large marginal sea between Greenland and Canada that formed during the early phase (Paleocene-Eocene) of North Atlantic opening (Gregersen et al., 2013;Oakey and Chalmers, 2012). Melville Bay is located in northeast Baffin Bay and overlies large areas of the northwest Greenland rifted continental margin (Fig. 1). The shelf area of Melville Bay is characterised by deep sedimentary basins, with thicknesses >10 km, separated by extensive elongate ridges Whittaker et al., 1997), with this complex paleo-rift topography overlain by a thick (1-3 km) Cenozoic post-rift succession ( Fig. 2) (Gregersen et al., 2013;Knutz et al., 2012Knutz et al., , 2015Knutz et al., , 2019.

Stratigraphic Framework
The rift and post-rift succession of the northwest Greenland continental margin has been mapped extensively and subdivided into seven seismic mega-units (mu) (A-G) bounded by regional horizons with updip unconformable expressions (Gregersen et al., 2013(Gregersen et al., , 2017Knutz et al., 2015) (Fig. 2). This regional seismic-stratigraphic framework for West Greenland has been adopted as a basis for a detailed 3D seismic analysis of the uppermost (1-2 km) stratigraphic packages.
The oldest section of the post-rift succession is represented by mu-E and consists mainly of hemipelagic marine mudstones intercalated with submarine fan deposits (Knutz et al., n.d). The mu-E succession was deposited during the seafloor spreading phase of Baffin Bay (syn-drift), followed by mu-D, of likely late Eocene to early-middle Miocene age. The lower part of mu-D consists of asymmetric wedges of submarine fans that possibly formed as a result of inversion as Greenland and North America converged and spreading in Baffin Bay ceased (Knutz et al., 2012;submitted). The upper part of mu-D displays a hemipelagic infill character, presumably with a high clay content facilitating the development of polygonal fault networks (Goulty and Swarbrick, 2005). Mu-D thins significantly above the crest of the Melville Bay Ridge (MBR) and is top bounded by a regional unconformity (horizon d1) (Fig. 2) (Cox et al., 2020a;Gregersen et al., 2019;Knutz et al., 2012). This thinning package represents the seal for a gas-charged mass transport deposit covering 420 km 2 of the MBR crest (within mu-D2) (Fig. 2). The reservoir displays seismic characteristics interpreted as direct hydrocarbon indicators and was likely deposited as a shallow marine spit complex along the ridge axis during the Eocene (Cox et al., 2020a). Reservoir deformation and transport occurred along a muddy décollement layer that represents horizon e1; an important unconformity that marks the top of the MBR structure and the transition from syn-to post-rift sedimentation .
On the inner shelf, the uppermost post-rift succession consists of extensive, thick marine sediment packages that were deposited during the late Miocene and Pliocene (mu-C and mu-B) (Knutz et al., 2015) ( Fig. 2), and have been heavily truncated by glacial erosion. Mu-C, which buried the MBR, is interpreted as a deep shelf drift accumulation influenced by northward flowing ocean currents (Knutz et al., 2015. Mu-B is considered mainly Pliocene in age and displays mounded and wavy depositional features that preferentially infilled erosional scarps where large parts of mu-C had been removed (Fig. 2). Basinward, the contourite accumulations and their correlative downslope mass transport deposits are buried progressively deeper by the thick glacigenic packages of mu-A corresponding to the Melville Bay Trough Mouth Fan   (Fig. 2). The progradational succession of mu-A has been subdivided into 11 units reflecting deposition by migrating ice-streams since the late Pliocene Newton et al., 2019).

Fig. 2.
Regional stratigraphy. (A) Uninterpreted regional 2D seismic reflection line in two-way-time, displaying the stratigraphy across the Melville Bay continental shelf. (B) An interpreted version of (A) defining the character of the regional stratigraphy including deep rift basins, the prominent Melville Bay Ridge and the thick package of glacigenic sediments. Nine seismic mega-units (mu-H to -A) as well as top unit horizons (horizon hx to b1) are defined, based on Gregersen et al. (2013Gregersen et al. ( , 2017 and Knutz et al. (2015). The occurrence of potential seismic fluid anomalies, faults and possible organic rich Cretaceous strata are also shown. The location of the seismic line is shown on Fig. 1.

Structural Framework
The paleo-rift topography consists of northwest-southeast trending elongate ridges formed during Early Cretaceous to early Paleogene synrift phases (Figs. 2 and 3) (Gregersen et al., 2013(Gregersen et al., , 2017Larsen et al., 2009;Whittaker et al., 1997). The Kivioq Ridge, located beneath the present outer shelf margin, and the MBR, separate the main depocentres of the Melville Bay Graben and South Kivioq Basin (Figs. 2 and 3) (Cox et al., 2020a;Gregersen et al., 2013Gregersen et al., , 2017Knutz et al., 2015). The ~200 km long MBR has continued as a positive relief feature into the Neogene, long after rifting ceased in the late Eocene/early Oligocene (Oakey and Chalmers, 2012;Welford et al., 2018). Consequently, the flanks of the MBR are onlapped and in parts draped by both syn-and post-rift sedimentary successions ( Fig. 2) with potential implications for regional fluid migration.
The MBR and Kivioq Ridge are characterised by numerous internal and ridge bounding deep tectonic normal faults generally oriented parallel to the ridge strike (Figs. 2 and 3) . The MBR formed as an elongate tilted fault block, with large offsets and disconnected syn-rift stratigraphy on its western flank adjoining the 'western ridge fault ' (Figs. 2 and 3). This is contrasted by a steeply dipping stratigraphy to the east that continues into the Melville Bay Graben (Figs. 2 and 3).
The deep pre-and syn-rift tectonic faults extend vertically to the horizon e1 unconformity marking the MBR apex ( Fig. 2)  . This horizon acts as a boundary between the deeper faults and a dense fault network within the post-rift stratigraphy (above horizon d2; Fig. 2). As well as tectonic faults, the 3D seismic mapping revealed that a large proportion of post-rift faults, especially within the grabens, exist as closely spaced, short offset faults that likely represent an intra-formational polygonal fault network formed by compactional processes, e.g. shale dewatering and diagenesis, (Figs. 2 and 3) (Cartwright and Dewhurst, 1998;Cox et al., 2020a). A neotectonic influence on the ridge complex is suggested by the observation that faults Fig. 3. Structural framework across the northern section of the study area. A two-way-time structure map of horizon e1 has been used to define the location of the Melville Bay Ridge. The location of deep tectonic faults (black) (including the western ridge fault) and post-rift faults (white) are shown, with the faults observed to strike parallel to ridge strike. Two windows (to scale) displaying seismic variance time slices from depths of 1600 ms TWT (Pitu) and 2230 ms TWT (Anu), display a network of thin, curvilinear high amplitude features that represent the polygonal fault network that exists within mu-D and -C, with the location of these faults across the area shown by the diagonal hatch. The location of Figs. 4, 9A, C, 15 and a section of 2 are shown.
intersecting mu-D and -C over the ridge flanks often extend from deeper structural features (Fig. 2).

Seismic and Velocity Data
The main database for this study consists of two 3D seismic reflection datasets and one regional 2D seismic reflection grid that covers much of the continental shelf in Melville Bay (Fig. 1). The first 3D seismic reflection survey (Pitu) was acquired by Cairn Energy in 2011 and covers an area of 1672 km 2 . It is provided in two-way time (TWT) down to 6.5 s and was acquired using 93 3D acquisition lines with 25 × 12.5 m bin spacing and a line separation of 1 km. A 25 m shot point interval (flip flop) was used with a 50 m source separation creating a common-midpoint fold of 70. Ten 7050 m long streamers were used, each with 564 channels and towed at a depth of 20 m. The dominant frequency varies with depth but is 28 Hz at approximately 1500 ms TWT, producing a dominant wavelength of 71 m (using an average velocity of 2.0 km/s) and a vertical resolution of 18 m. The reprocessing of this survey in 2013 provided a sub-set high resolution 3D seismic survey in TWT (Pitu HR survey) that covers an area of 1135.5 km 2 down to 5 s TWT. The resolution of the data has been improved to a spatial resolution of 12.5 m (inline) x 6.25 m (crossline) and a vertical resolution of 11 m at the same stratigraphic interval (~1500 ms TWT) due to a dominant frequency of 45 Hz and a dominant wavelength of 44 m.
The processing of the Pitu survey involved noise attenuation and filtering, 3D deconvolution, Pre-Stack Depth Migration (PSDM), postmigration velocity analysis and a parabolic radon de-multiple process. An understanding of the velocity field was gained during the PSDM of the full survey, which involved an iterative velocity investigation including multi-layer tomography to minimise residual move-out, stack scanning to help flatten irregular unconformities and finally Kirchhoff pre-stack Tilted Transverse Isotropy (TTI) depth migration, creating a depth domain volume which was subsequently converted into the time domain. This process produced an interval velocity cube based solely on seismic velocities (Fig. 4) as well as a depth domain seismic volume down to 10 km. An assessment of the accuracy of both the depth conversion and interval velocity cube is not possible as no calibration data (such as boreholes) exist in the area for comparison. However, interval velocity variations show a good correlation to both regional 2D stacking velocities and the location of regional unconformities (where rapid increase in velocity may be expected), such as at the top of the MBR (horizon e1) (Fig. 4).
The second 3D seismic reflection survey (Anu) was acquired by Shell in 2013 over an area of 8700 km 2 (Fig. 1). It is provided in TWT down to 7.5 s and was acquired using 118 sail lines with a 50 m (cross-line) x 6.25 m (inline) bin spacing and a line separation of 600 m. The acquisition used a dual vessel set-up with six, 7050 m long streamers (each with 564 channels) that were separated by 200 m and towed at a depth of 10 (front end) to 15 m (tail end). Two sources were separated by 100 m at a depth of 8 m and were fired using a shot point interval (flip flop) of 25 m, producing a nominal fold coverage of 70. The processing workflow included noise and acquisition footprint attenuation processes, multiple removal using both 3D surface related multiple reflection (SRME) and 2D model-based water-layer de-multiple (MWD) modelling, a 1 km migration velocity analysis and a final isotropic Kirchhoff pre-stack time migration (PSTM). This produced a volume with a dominant frequency of 30 Hz at approximately 1500 ms TWT, a dominant wavelength of 67 m (using an average velocity of 2.0 km/s) and a vertical resolution of 17 m.
The 15 lines of 2D seismic reflection data used were acquired by TGS between 2007 and 2009 and form part of the Baffin Bay 2D (BBRE11) regional dataset (Fig. 1). The lines used here cover 5038 km and were acquired using a single 6000 m streamer with 480 channels, 12.5 m intervals and a sampling rate of 2 ms, as well as a 2000 psi source with a 25 m interval at a gun depth of 8 m, yielding a 120-fold stack. The survey was reprocessed in 2011 to boost low frequencies and enhance resolution, specifically targeting multiple and bubble pulse attenuation. The processing sequence involved a 2 km Kirchhoff pre-stack curved ray time migration and velocity analysis within this process produced stacking velocity data for each 2D line.
All data used were provided in SEG normal polarity with a downward increase in acoustic impedance (such as at the seabed) represented by a red positive peak and a downward decrease in acoustic impedance represented by a blue negative trough (Fig. 5). No well data were available within the area of interest.

Interpretation Methods
The regional structural and stratigraphic framework (Gregersen et al., 2013(Gregersen et al., , 2017Knutz et al., 2015) was propagated through the 3D seismic data using Schlumberger's Petrel software and standard 3D seismic interpretation techniques (Cox et al., 2020b;Posamentier, 2004;Posamentier et al., 2007). This study focused on the mapping of features such as seismic discontinuities and amplitude anomalies that could be linked with fluid migration and hydrocarbon accumulations (Hilterman, 2001). Analysis of such features found that likely gas-related anomalies Fig. 4. Seismic cross section from the Pitu HR survey demonstrating the relationship between seismic stratigraphy, interval velocity and potential gas hydrate deposits. The seismic amplitudes are overlain by a semi-transparent cross section displaying the corresponding interval velocity. Key seismic horizons and a bottom simulating reflector (BSR), marking the base of the gas hydrate stability zone (GHSZ), are indicated. 1D interval velocity profiles display the variation of velocity through the GHSZ and the adjacent stratigraphy. The cross section location is shown on Fig. 3. often displayed seismic amplitudes larger than − 10,000 (Fig. 5), although dimmer anomalies may also comprise shallow gas accumulations associated with saturations of less than 10%, poorer reservoir quality or thin-bed effects (Hilterman, 2001), and this was considered when filtering results. Using this understanding of anomaly amplitude, automated anomaly detection and filtering was applied to the 3D data via the method described in Cox et al. (2020c) (Fig. 5). The individual scrutiny of each anomaly within the filtered results was then conducted to determine whether that anomaly was in fact fluid-related or created by another feature.
Non-fluid-related bright 'soft' amplitude anomalies may represent features such as high porosity, low density and low velocity (i.e. low acoustic impedance) stratigraphic layers such as diatomaceous ooze or coal beds, organic-rich claystones, and data acquisition footprint (stronger near the seabed) (Cox et al., 2020c). Distinguishing between fluid-related and lithology-related anomalies, for example those created by low density oozes which may create similar amplitudes, is sometimes difficult, but can usually be distinguished based on context as lithological variations are more often stratigraphically restricted and more widespread, e.g. Batchelor et al. (2017). This evaluation is therefore somewhat subjective, increasing uncertainty, and this should be considered in the final result. Additionally, high saturation gas anomalies have likely been detected efficiently through this process, but the number of fluid-related anomalies may be underestimated, as very thin (few metres) or very low saturation gas pockets (<<5%), causing only low amplitude seismic anomalies, may have been omitted.

Geothermal Gradient and Heat Flow Calculations
Near surface geothermal gradient and heat flow were estimated across part of the Pitu 3D seismic survey area using an identified BSR that is thought to represent the base of the gas hydrate stability zone (GHSZ) (Dickens and Quinby-Hunt, 1994;Grevemeyer and Villinger, 2001). Gas hydrate stability is influenced by fluid composition, pressure and temperature conditions (Kvenvolden, 1993), and the estimation requires key parameters including temperature at the seabed and BSR, as well as the thermal conductivity of the shallow gas hydrate hosting sediments (Grevemeyer and Villinger, 2001;Minshull and Keddie, 2010). These parameters are often provided by seabed and downhole temperature probes and logging data. However, no such data exist in northeast Baffin Bay, so these essential parameters have to be either extracted from interpolated global databases or from published gas hydrate stability models and empirical relationships using the constraints provided by the seismic data (TWT and depth of seabed and BSR, interval velocity of the GHSZ).
Estimating parameters in such ways creates uncertainty -e.g. estimating heat flow without direct temperature data is thought to cause an error of 20% or above, whilst the estimation of thermal conductivity through empirical relationships instead of seabed measurements, can cause errors of 5-30% (Grevemeyer and Villinger, 2001;Minshull, 2011). Other uncertainties arise from seismic interpretation errors, velocity model creation and depth conversion and uncertainty of the true phase boundary depth, due to an unknown composition of pore water and hydrate (other gases such as hydrogen sulphide or higher order hydrocarbons can affect the phase boundary) (Grevemeyer and Villinger, 2001;Minshull, 2011;Minshull and Keddie, 2010). Error resulting from these uncertainties is often over 10%, with the compound uncertainty being potentially large (e.g. > 50% cf. Grevemeyer and Villinger, 2001), but in areas where no data exist, an estimate with recognised uncertainty is better than no information at all. The geothermal gradient and heat flow were estimated using the procedure described by several authors (Dickens and Quinby-Hunt, 1994;Grevemeyer and Villinger, 2001;Yamano et al., 1982) and used within a number of similar studies (Calves et al., 2010;Minshull and Keddie, 2010;Serié et al., 2017;Shankar et al., 2010).
It was assumed that hydrate gas is pure methane and the pore water exhibits normal sea-water salinitythis is the most commonly observed scenario and thus a typical assumption when no borehole data are available (Dickens and Quinby-Hunt, 1994;Shankar et al., 2010). There is a chance of higher order hydrocarbons, with the hydrate forming gas potentially sourced from Cretaceous marine source rocks (Cox et al., 2020a;Nohr-Hansen et al., 2018), however, no 'double-BSR' is observed that may suggest the presence of three fluid phases Geletti and Busetti, 2011). The first step was to interpret the seabed and the BSR on the Pitu depth domain seismic. The pressure (ρ) at the BSR was then calculated using Eq. 1, which represents a standard hydrostatic pressure equation (Dickens and Quinby-Hunt, 1994). This used an estimated seawater density of 1034 kg/m 3 for relatively deep,

Fig. 5.
A seismic cross section (A) illustrating the method of using proportional minimum amplitude extraction windows to identify shallow gas anomalies. Examples of isolated anomalies are shown (category I). The polarity of seismic data used in this study is shown with an increase in acoustic impedance at the seabed representing a positive amplitude brown-red-yellow peak. The location of the line is shown on Fig. 7. (B) Histogram of the amplitude distribution from an extracted minimum amplitude window, demonstrating the filtering process for isolating the seismic signal most likely representing shallow gas occurrences. saline and cold ocean bottom waters on the Melville Bay shelf provided by Tang et al. (2004). The pressure equation assumes hydrostatic pressure at the BSR depth and therefore depth is measured from mean sea level.
BSR Pressure (ρ) = (1034*9.81*BSR Depth)/ − 1000000 (1) Units: Seawater density = 1034 kg/m 3 Gravity = 9.81 m/ s 2 Megapascals (MPa) conversion = 1,000,000 The temperature at the BSR was estimated from the calculated pressure using Eq. 2. This equation is based on an empirical relationship between pressure and dissociation temperature, defined by experimental data for various methane hydrate stability conditions, with a regression showing a Chi-squared of >0.99 (Dickens and Quinby-Hunt, 1994). Eq. 2 represents this empirical relationship rearranged to find the dissociation temperature for any pressure between 2.5 and 10 MPa. The only variable in the equation is pressure, with the values 0.00379 and 0.000283 representing constants derived from the original empirical relationship.
BSR Temperature = (1/((0.00379)-(0.000283*Log(BSR Pressure) ) ) )-273 (2) The temperature change across the GHSZ was then calculated by subtracting the seabed temperature which was estimated from the world ocean temperature database (Locarnini et al., 2013), showing an annual mean temperature of ~1.5 • C for ocean bottom waters at this depth (~650 m) across Melville Bay. The seabed elevation varies in depth by ~100 m across the BSR area and therefore, this temperature may vary slightly (less than 0.5 • C), translating to a BSR temperature error of < +/− 0.2 • C. The temperature difference across the GHSZ was then divided by the thickness (in metres) (Eq. 3) to calculate the Geothermal Gradient (Dickens and Quinby-Hunt, 1994;Grevemeyer and Villinger, 2001).
Units: BSR/Seabed Temperature = • C BSR/Seabed Depth = m Heat flow across the same zone was then derived from the estimated geothermal gradient. This requires an understanding of the thermal conductivity structure of the GHSZ. In the absence of borehole control, an empirical relationship between interval velocities and thermal conductivity (k) is used, that was derived from a number of experimental datasets (filtered to only include wet samples to represent subsurface fluid saturation) that include P-wave velocity, thermal conductivity and porosity measurements across a wide range of lithologies, e.g. Grevemeyer and Villinger (2001), Boulanouar et al. (2013) and Esteban et al. (2015). This relationship was used to estimate the thermal conductivity of the shallow sediments as a function of the seismically-defined interval velocities of the GHSZ (Eq. 4). The empirical relationship is represented by 0.5071 in eq. 4, which is the intercept of the regression equation through the selected experimental data points. Within eq. 4, the interval velocities are multiplied by 0.001 to convert from m/s to km/s to ensure the thermal conductivity units match those of geothermal gradient and heat flow.
Thermal Conductivity (k) = (Interval Velocity*0.001)-0.5071 Units: k = W/m/K Interval Velocity = m/s The final step was to use the thermal conductivity and the geothermal gradient to calculate heat flow (Q) through the GHSZ using Fourier's Law (Eq. 5). It should be noted, as acknowledged above, that the compound uncertainties may be as much as 50-60% on the final heat flow estimates (Grevemeyer and Villinger, 2001).

Shallow Gas
Numerous, scattered seismic anomalies are observed within the Cenozoic sedimentary succession covering the rift basins of Melville Bay (Figs. 6 and 7) and potentially represent the presence of free gas or gasrich pore fluids. The identified acoustic anomalies often exhibit a bright, negative amplitude top reflection against a generally lower amplitude background response (Fig. 5). The anomalies vary in their seismic character and have been subdivided into four categories:

Category I: Isolated anomalies
The majority of seismic anomalies observed across the study area occur as isolated, negative amplitude anomalies within a background of subtle or low amplitudes (Figs. 6, 7 and 8). These isolated (category I) anomalies still vary in seismic character and therefore have been subdivided into four key types (Figs. 5 and 7B-D; Type 1-4). The vast majority of category I anomalies display a small, bright event, characterised by a single loop trough-peak response (Type 1) (Fig. 5). These anomalies are limited in the horizontal (<~150 m) and vertical (<~20 m) plane but tend to occur at various stratigraphic levels (Figs. 5, 7B and 8). Polarity reversals are occasionally observed at the edge of the anomaly, dependent on the seismic properties of the host sediments.
Isolated anomalies also occur as larger features observed across multiple stratigraphic layers (horizontal/vertical extent of >2 km/ >200 m) characterised by bright negative amplitude top reflections and a positive amplitude base (Type 2) (Figs. 5 and 7A). Internal reflections exhibiting bright amplitudes, as well as polarity reversals at the edge of the stacked anomalies, suggests a variable pore fluid as the main cause for the lower seismic velocities (Fig. 5).
Within mu-D1, above the southern extension of the MBR, a concentration of category I anomalies occurs (Type 3). These anomalies differ from the more typical Type 1, as they appear restricted to a single stratigraphic horizon, with the anomalies and the surrounding stratigraphy being offset by closely spaced, near vertical linear features characterised by dim amplitudes (Figs. 7 and 8A). These offsets are interpreted as small faults that extend downwards to the base of the Eocene gas reservoir (see category IV). A seismic variance attribute extraction through the stratigraphy at this depth shows a dense network of these faults ( Fig. 9A-B). The bright anomalies occur in between the faults suggesting entrapment of hydrocarbon fluids focussed within the Miocene interval directly above the MBR (Figs. 8 and 9A-B).
Several, more laterally continuous faults are also observed on the variance extraction (Fig. 9), with some containing small, bright amplitude seismic anomalies within or truncating against the low amplitude linear zone between the offset horizons, a zone likely representing the fault plane (Figs. 8 and 9D).
Numerous category I anomalies are observed at the base of mu-C in a limited area in the southern part of the Pitu HR survey (Fig. 7A purple polygons). These anomalies, referred to as Type 4, occur as vertically stacked, negative amplitude narrow anticlines across several layers of stratigraphy. The anticlinal anomalies follow the structure of the host stratigraphy, and often gradually diminishing upwards suggesting drape across an underlying structure (Figs. 7D and 8B). In the lower section of mu-C (Fig. 8B), the anomalies occur in between fault locations within the underlying and more chaotic section of mu-D. They also become more widespread (~2 km) above a local unit of enhanced thickness within mu-D1 (Fig. 8B), which causes the overlying anomalies and the surrounding stratigraphy to dip more steeply around its crest.

Category II: Anomalies terminating at the base GHSZ
This category refers to bright, negative amplitude anomalies observed within dipping strata packages that terminate at a level D.R. Cox et al. Fig. 6. A regional map of the area covered by 3D seismic data, displaying all seismic anomalies that were identified and interpreted to represent shallow hydrocarbon occurrences, the location of paleo-rift topography and the direction of glacial unit progradation. Location of Figs. 7 and 11 are shown. corresponding to a BSR (Figs. 7 and 10). The BSR appears as a crosscutting negative amplitude reflection, interpreted as the base of the GHSZ (see section 4.2) (Figs. 10). These anomalies exist within multiple layers of tilted stratigraphy with the bright negative response existing for ~50 ms below the BSR (Fig. 8C). The brightest terminating anomalies exist within linear zones that trend northwest-southeast (Figs. 9C and 10C) and represent the location where the base GHSZ cross-cuts the host layers. These host layers maintain a bright reflection down-dip from the truncated anomaly, but switch to a positive amplitude response that continues regionally across the study area (Fig. 10A). The category II anomalies are most clearly expressed above laterally and vertically extensive faults that extend from the top horizon (e1) of the western flank of the MBR, to a near-seabed position (Figs. 2, 3, 9C and 10).

Category III: Steeply dipping bright horizons
This anomaly category is observed within the progradational units of mu-A, forming part of the Melville Bay Trough Mouth Fan (Figs. 1, 2, 6 and 11) . It is characterised by anomalously bright, negative amplitude top reflections, that display top set erosion and are interpreted as glacial clinoform packages (Fig. 11) (Newton et al., 2020). The clinoforms covered by the Anu 3D seismic survey display west to south-westerly dips (Figs. 6 and 11). The category III anomalies are distributed throughout much of the shallow stratigraphy (uppermost 1 s TWT) within the Anu 3D survey (light green area; Fig. 6), but with the highest amplitudes observed below the truncated top-sets, representing paleo-shelf breaks of the prograding wedge (dark green areas; Fig. 6). Towards the northwest within the survey area, corresponding to the oldest prograding units, dipping bright horizons terminate abruptly updip against a major glacigenic unconformity (Fig. 11). Down-dip, the category III anomalies continue along semi-continuous, seismic-stratigraphic horizons but with variable and generally fading amplitudes (Fig.  11).

Category IV: Eocene gas reservoir
An extensive, potential gas reservoir has been mapped at the base of the Cenozoic succession (mu-D2) overlying the MBR (Fig. 8) (Cox et al., 2020a). The reservoir is likely of Eocene age and interpreted as laterally continuous sand deposits on the northern parts of the ridge crest while towards the south and east, the unit appears fragmented due to mass transport (Figs. 8 and 12). Thus, large sections of the reservoir exist as separated blocks above a likely muddy decollement surface (represented by horizon e1) (Figs. 2 and 8). The reservoir package displays very bright amplitudes on seismic data including a bright negative amplitude top reflection and a bright positive amplitude base (Fig. 8a).
The thin overlying sealing stratigraphy of mu-D1 is characterised by a dense network of closely spaced small offset faults that often appear aligned to the edge of the reservoir blocks (Figs. 8 and 9A-B). These faults, hosting isolated bright anomalies (category I), commonly extend between the base of the reservoir (horizon e1) and horizon d1, with some extending upwards towards the seabed. Hence, there is a clear spatial relationship between the regional distribution of the reservoir blocks on the MBR crest and the isolated anomalies throughout the shallow stratigraphy (mu-D1 and -C) (Figs. 6 and 12).

Gas Hydrates
A discontinuous BSR has been observed across an area of 537 km 2 above the central axis and western flank of the MBR (Figs. 6, 7 and 10). This feature exists approximately 200 m below the seafloor, in water depths of 625-720 m and marks the base of the GHSZ. The BSR is observed as a negative amplitude reflection with variable intensity that cross-cuts the stratigraphy. It is best defined in areas with tilted stratigraphy where bright anomalies terminate at the BSR and the base of the GHSZ (category II) (Figs. 7E and 10B). Away from these anomalies, and especially in areas where the stratigraphic dip is similar to the seabed the BSR amplitude is reduced to a dimmer response, but can still be observed (Figs. 8, 9D-E and 10). The seafloor topography across the region contains several erosional and depositional features related to glaciation (Figs. 1 and 14). The BSR is observed to mirror this seafloor topography, deepening and shallowing in areas of erosion and deposition, respectivelye.g. such as the glacial wedge (Figs. 10A and 13G-H). The thickness of the GHSZ varies, generally in relation to seabed depthe.g. with the GHSZ thinning or thickening in areas of shallower or deeper water, respectively. An exception to this trend is observed in the northern and eastern parts of the BSR, where the GHSZ thins in areas where water depth increases ( Fig. 13C and E). The increased water depths are a result of seabed glacial erosion (Fig. 14)  .
Interval velocities increase by ~200 m/s within the stratigraphy above the BSR (e.g. GHSZ), when compared with areas outside the BSR (Fig. 4). Below the BSR, seismic velocities at the GHSZ base are reduced from ~2200 m/s within the GHSZ to ~1800-1900 m/s beneath. Conversely, a velocity increase is observed at a corresponding depth level away from the BSR (Fig. 4). Despite its discontinuous appearance, the distribution of the BSR suggests a relationship to the paleo-rift topography and the MBR. The BSR is observed in most areas along the western flank of the MBR, but disappears above a deep valley between the northern and southern extensions of the ridge (Fig. 10A). In the vicinity of the valley system, the amplitude strength of the BSR reduces, and disappears completely above the central parts of the valley. Within the buried valley zone, several of  Fig. 3. D) A seismic cross section in two-way-time from the Pitu 3D survey displaying potential fluid anomalies trapped within and truncated against a vertically extensive tectonic fault. This fault is observed to cross cut a low amplitude section of the bottom simulating reflector (BSR), which is also shown on (E). The location of both (D) and (E) is shown on (C). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) the regional horizons that host the bright category II anomalies are located at depths well below the depth of the base GHSZ. The lack of cross-cutting between the host strata and the GHSZ within this area is a likely reducing factor for the seismic imaging of the BSR (Fig. 10A and C). In areas containing extensive faults that extend from the top ridge horizon (horizon e1) and through the post-rift stratigraphy, the BSR retains a strong amplitude response within the host strata (Figs. 2, 3, 9C and 10).

Evidence for Seabed Fluid Escape
The seafloor located in the northern part of the study area represents a major glacial trough, characterised by numerous features of glacial erosion and deposition, such as lineations, ridges, iceberg scours and grounding zone wedges (GZW) (Figs. 2, 13 and 14) . Within the trough area, semi-circular seafloor depressions are commonly observed, with lateral dimensions up to 500 m and depths up to 50 m (Fig. 14D-B). Small mounds of sediment often occur on the seabed next to the depression. These mounds are not evenly distributed but tend to be skewed to one side of the depression (Fig. 14D).
The depressions seem to be unrelated to the underlying stratigraphy, as no seismic features are observed consistently beneath their locations, such as faults, pipes or anomalous brightening (although on occasion a fault will occur beneath the depression, e.g. Fig. 14E). On the seabed structure map, curvilinear scours are often observed connecting the circular depressions, commonly extending towards the northeast, trending similar to glacial lineations ( Fig. 14B and D). In the western part of the study area (Anu survey), numerous seabed depressions occur towards the southwest, lee-side of a depositional feature interpreted as a GZW . Several semi-circular depressions are also observed on the base horizon of the GZW, which represents a paleoseabed surface (white dashed line on Fig. 14C). Both the seabed depressions and GZW deposits coincide with lower seismic amplitudes of the seabed horizon (Fig. 14E). Based on seabed geomorphology and the distribution of glacial features, the semi-circular depressions are interpreted as having been formed by glacial erosion, and most likely A discontinuous BSR can be observed above much of the paleo-rift structure, but is not observed above the majority of the valley area. The BSR varies in amplitude between areas containing terminating bright anomalies at the base gas hydrate stability zone (GHSZ) and areas without these anomalies where it displays as a low amplitude reflection (highlighted by white brackets). B) A seismic cross section in two-way-time from the Pitu HR 3D survey showing a clear cross-cutting BSR, as well as numerous bright anomalies terminating at the base of the GHSZ (category II). The location of both (A) and (B) are shown on (C). C) A minimum amplitude extraction for a 20 ms two-way-time window across the mapped BSR horizon that highlights the location of bright amplitude anomalies that terminate at the base of the GHSZ (category II). The relationship of the hydrate deposits to the underlying Melville Bay Ridge is also shown, via the underlying two-way-time structure map of horizon e1 defining the ridge topography. The intersection of a regional potentially porous horizon (identified on A) with the base GHSZ is also shown (yellow dashed line), with the intersection location coinciding with numerous terminating fluid anomalies. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) represent iceberg pit-marks (e.g. Brown et al., 2017) as opposed to gashydrate explosion craters (e.g. Andreassen et al., 2017).

Geothermal Gradient
The geothermal gradient has been estimated using the GHSZ thickness inferred from the BSR distribution covered by the Pitu 3D survey (depth domain) and known P-T phase relationships for gas hydrates ( Fig. 13; section 3.3). Accordingly, pressure and temperature at the BSR depth, range between 8.0 and 9.5 MPa (average 8.9), and 10.0-11.6 • C (average 10.9) respectively (Fig. 13D). As both the pressure and temperature estimations are dependent on BSR depth, the two parameters show identical lateral variation. Hence, shallowing of the BSR coinciding with a shallowing seabed depth, as seen in the southern part of the Pitu survey (in the area of glacigenic deposition), causes a reduction in both pressure and temperature (Fig. 13). A differing trend to that of the seabed depth is seen in the northern section and along the eastern edge where the seabed deepens, at the same locations where the relationship between seabed depth and GHSZ thickness diverts (section 4.2).
The calculated geothermal gradient across the thickness of the GHSZ ranges from 40 to 59 • C/km with an average of 49 • C/km (Fig. 13). The trend of the estimated geothermal gradient matches that of the GHSZ thickness, with the geothermal gradient increasing in areas of a thinner stability zone and vice versa. However, superimposed on this trend, are a number of localised high values that are either concentric or linear in shape. These anomalous values occur in the same location as similar features observed on both the GHSZ thickness (Fig. 13E) and the seabed depth maps (Fig. 13C).
The interval velocities through the GHSZ range between 1950 and 2350 m/s (Fig. 13F) and are observed to increase towards the south within an area coinciding with glacigenic deposition at the seabed (Fig. 14 C, and F -H). The conversion of these velocities to thermal conductivity yields a range from 1.45 to 1.85 W/m − 1 /K − 1 . The thermal conductivity estimates, in combination with the geothermal gradient model values, produces a heat flow across the thickness of the GHSZ, ranging from 58 to 100 mW/m 2 with an average of 81 mW/m 2 (Fig. 13B). Again, several concentric or linear anomalies are observed on the heat flow map (Fig. 13B), which coincide with areas of high geothermal gradient, a thin GHSZ and high interval velocities.

Hydrocarbon Occurrences
The detailed analysis of the shallow Cenozoic stratigraphy in Melville Bay yielded numerous seismic amplitude anomalies which have been interpreted to represent the most northerly occurrence of shallow gas and gas hydrates offshore western Greenland (Figs, 6-8 and 11). The gas anomalies were categorized by seismic character and mainly represent the different trapping mechanisms observed across the region (categories I-IV), with the majority of hydrocarbon occurrences existing as isolated pockets of gas (category I), with variable seismic character (Type 1-4) (Figs. 5-9).
Free gas was observed trapped against the base of the GHSZ (category II) (Figs. 6-10), a boundary recognised by a discontinuous BSR (Berndt et al., 2004;Hillman et al., 2017), hosted within tilted, likely porous horizons, with the termination of the fluid response down dip (Fig. 10B), suggesting a maximum gas column of ~50 m. Dimming of the BSR (Figs. 8, 9D, E and 10) occurs in areas away from these terminating free gas anomalies, as the acoustic impedance contrast at the base GHSZ is lowered (Hilterman, 2001). The free gas may be absent due to a lack of migration, lack of reservoir at the base GHSZ, and due to all available gas having been converted to hydrate. The dimming may suggest low saturation of gas hydrate in areas without a BSR (Spence et al., 1995). It is also possible that dim or absent BSRs could be a result of recent glaciations having caused variations in the temperature field which still persist today, if the hydrate system has not fully re-equilibrated postglaciation (Crémière et al., 2016;Grassmann et al., 2010;Mienert et al., 2000). Overall, the abrupt alternations between BSR and no BSR seen in Fig. 10 are probably due to a combination of stratigraphic occurrence of reservoir porosity and permeability and availability of methane around the base of the GHSZ.
Furthermore, the gas interpreted to exist within glacigenic, progradational units of mu-A (category III) , is observed trapped up-dip against horizontal and unconformable muddy strata of the overlying glacial sequence (Fig. 11). The brightest fluid-related amplitudes coincide with truncated top-sets at the paleo-shelf break positions of the prograding wedge; possibly representing increased gas concentrations or greater porosities (Hilterman, 2001). The rotating alignment of these fluid anomalies southward, as well as the variable dip of glacigenic stratigraphy (Fig. 6), reflects the variable pathways of ice streams and shelf break progradation Newton et al., 2019Newton et al., , 2017.

Fluid Migration
Although the seismic evidence provided only offers a static view of contemporary hydrocarbon occurrences, and features such as chimneys, pipes and pockmarks are not observed at present day (Chand et al., 2012;Huuse et al., 2010), several relationships between structural and stratigraphic elements and the distribution and character of the hydrocarbons observed, provides evidence towards understanding the migration history of fluids within this area.

Paleo-rift topography and tectonic faults
A significant concentration of free gas anomalies, as well as all identified gas hydrate deposits, exist directly above or along the western flank of the underlying MBR (9)(10), with the deposits only becoming discontinuous above a large buried valley separating the ridge's north and south extensions (Figs. 6 and 10). This relationship suggests that the underlying paleo-rift topography is influencing upward fluid migration. Post-rift depositional units either onlap onto the ridge flanks, or extend and thin above its crest, forming a broad anticline structure throughout mu-D and -C (Fig. 2) Gregersen et al., 2013). This causes the post-rift stratigraphy within the adjacent grabens to tilt towards the ridge structure, likely facilitating the up-dip migration of fluids into areas above the MBR.
Potential Cretaceous source rock intervals within mu-F (Nohr-Hansen et al., 2018; Planke et al., 2009), are interpreted to extend from the Melville Bay Graben up-dip into the MBR structure on the eastern flank; a factor which likely led to the charging of the Eocene reservoir (category IV) (Cox et al., 2020a). On the western flank the syn-rift stratigraphy terminates against deep tectonic faults (Figs. 2, 3 and 15), with disconnection of this strata occurring at the western ridge fault. This termination and offset likely restricts fluid migration from the South Kivioq Basin, and instead, fluids are forced upward through a number of post-rift fault conduits that are observed to extend from the tips of deeper faults and the top ridge horizon (horizon e1), to close to the seabed (Fig. 2-3, 9 and 15). Above these faults on the western ridge flank, gas hydrate deposits (Figs. 6,(9)(10) suggest the location of these migration pathways have influenced gas hydrate formation.
Furthermore, free gas trapped beneath the hydrate deposits (category II) is also observed above these fault locations (Fig. 9C), although this distribution may instead be influenced by the intersection of the host strata (regional, potentially porous and permeable horizons), that have been uplifted and tilted by the MBR structure, with the base GHSZ along the western flank of the ridge (Figs. 9C, 10 and 15). Still, no free gas is observed within the GHSZ above the fault locations (Figs. 10 and 15) and no seabed expulsion features are observed directly above the fault tips (Figs. 10 and 14-15). Therefore, any fluids migrating upward through the faults, post hydrate formation, either became trapped within the fault plane at the base GHSZ (Fig. 15), or migrated stratigraphically up-dip once reaching the potentially porous carrier bed horizons until again becoming trapped at the base GHSZ (Figs. 10 and 15).

Fig. 12.
A two-way-time structure map of horizon e1 that represents the top of the Melville Bay Ridge overlaid by the distribution of the Eocene gas reservoir (category IV) across the ridge crest (black polygons) (Cox et al., 2020a(Cox et al., , 2020b. The location of Isolated anomalies (category I) (red polygons), interpreted to represent shallow gas, are shown with their distribution following the trend of the underlying reservoir. This spatial relationship may suggest that the deeper reservoir is leaking hydrocarbons into the overlying stratigraphy. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Reservoir leakage
Above the crest of the MBR, a clear relationship is observed between the location of isolated gas anomalies (category I) and the trend of the Eocene gas reservoir (category IV) (Figs. 6, 8 and 12). The sealing unit (mu-D1) is at its thinnest above the reservoir (due to depositional thinning above the ridge and erosion) (Cox et al., 2020a;Gregersen et al., 2013), and is observed to contain a dense network of short offset, near vertical faults, extending from the base of the reservoir (horizon e1) to the top of the sealing unit (horizon d1) (Cox et al., 2020a). Several fluid anomalies are observed trapped within this dense network of faults (Figs. 8A and 9A-B), with their occurrence, as well as overlying anomalies in the lower section of mu-C (Figs. 8A, 10A and 12), suggesting that the underlying reservoir is leaking hydrocarbons upwards through these fault conduits (Fig. 8-9, 12 and 15) Ingram and Urai, 1999).
These faults occur in a dense, polygonal network (Figs. 3 and 8-9) and are geometrically similar to sets of polygonal faults within the same stratigraphy above the grabens (Figs. 2 and 15) (Cartwright and Dewhurst, 1998;Cox et al., 2020a). However, their exact origin is uncertain and fault growth in certain areas may instead be controlled by the presence of the underlying deformed reservoir (Figs. 8B and 12), with normal faults forming along the flanks of reservoir blocks, likely due to differential compaction above the blocks and adjacent areas (Fig. 8) (Cox et al., 2020a). Polygonal faults are then thought to nucleate at the tip of these faults (at horizon d1) and extend vertically into mu-C. A similar scenario may have occurred above the ridge crest, where the gas anomalies are observed (Figs. 8-9). Here, the underlying reservoir is more laterally continuous, but faults may still be influenced by reservoir structure as they coincide with irregular offsets between narrow sections of the reservoir (Fig. 8A), although this offset may instead be caused by the faults themselves. Therefore, if fault genesis is not typically polygonal (Cartwright et al., 2003), then gas leakage through these faults is more likely as polygonally faulted muds are known to create excellent sealing formations (Cartwright, 2019;Goulty, 2008). If leaking did occur, then this may be providing at least part of the methane constituent of the overlying hydrate deposits (Figs. 7-10).

Glacial clinoforms
The source and likely migration pathway of the gas interpreted to exist within glacial clinoform units of mu-A (category III) (Figs. 2 and   Fig. 13. Maps showing the calculated average geothermal gradient (A) and average heat flow (B) through the gas hydrate stability zone (GHSZ) interval for the area containing the observed bottom simulating reflector (BSR) within the Pitu depth survey limits. Maps used within this calculation are shown including seabed depth (C), BSR depth, pressure and temperature (D), GHSZ thickness (E) and average interval velocity through the GHSZ (F). Seismic cross sections in two-way-time (G) and depth (H) are shown across a section within the south of the BSR distribution highlighting the irregular seabed topography, and the likely post-glacial re-equilibration of the hydrate phase boundary beneath a depositional glacial wedge. 11) is uncertain, due to the disconnection of this stratigraphy from potential source rock horizons in the deeper basins ( Fig. 2) (Cox et al., 2020a;Nohr-Hansen et al., 2018) via thick packages of likely sealing formations including hemipelagic muds (mu-D) and contourites (mu-C and -B) (Gregersen et al., 2013;Knutz et al., 2015). No deep tectonic fault connections are observed, and the paleo-rift topography close to the shelf edge ( Figs. 1 and 2) is likely too far away to focus fluids into sediments that were only deposited since ~2.7 Ma . Shallow source rocks may exist, although sufficient burial and time for maturation is unlikely when considering burial only increased after ~2.7 Ma Wilson, 1975). Instead, the gas could be biogenically generated, either in situ as the clinoforms may contain organic material -as suggested for other glacial related clinoforms, e.g. Verweij et al. (2012);Muller et al. (2018) or more likely, from organic horizons within the contourite succession, e.g. Knapp et al. (2019); Rebesco et al. (2014). The contourites (mu-C and -B) are estimated to exist at temperatures around 40-50 • C at present day (Fig. 11), the theoretical peak temperature of biogenic generation (Stopler et al., 2014), based on the estimated shallow geothermal gradient from this study (49 • C/km) (Fig. 13).

Seabed expulsion
Attempts were made to document fluid expulsion features at the seabed (e.g. chimneys, pipes and pockmarks) which may represent the existence of cold hydrocarbon seeps into the water column Chand et al., 2012), similar to those observed ~500 km to the south near Disko Island (Nielsen et al., 2014). The seafloor across the study area, however, represents a time-transgressive glacial unconformity that has experienced significant erosion since ~2.7 Ma (Figs. 2 and 14) (Bennett et al., 2014;Knutz et al., 2019;Newton et al., 2019). This erosion has removed large amounts of sediment from the inner shelf and likely removed any seabed expulsion features that may have formed prior to glaciation (Fig. 14). Further towards the shelf edge, horizons interpreted as pre-glacial paleo-seabeds are now buried beneath a thick unit of glacigenic sedimentation (Fig. 2), and seismic resolution is likely insufficient to distinctly image any potential buried features. Additionally, much of the area above the MBR, where upward fluid migration is focussed, contains gas hydrates (Figs. 2,(6)(7)10 and 15), which may restrict upward fluid migration and expulsion to the surface (through both conventional trapping -such as the free gas anomalies observedand transformation to hydrate) (Grauls, 2001;White, 1979).
Several erosional depressions are however, observed on the present day seabed within two distinct regions of the 3D data (Fig. 14). These features could represent fluid escape locations at the seabed Chand et al., 2012), but no direct evidence for fluid flow is observed, although features such as gas chimneys would likely have had time to dissipate if formed during the early stages of glaciation . Evidence such as their geometry, the asymmetrical sediment berms and their occurrence within seabed scours (Fig. 14), instead suggests that the depressions are more likely associated with glacial erosion of the seabed, and are interpreted as iceberg pit-marks (Brown et al., 2017). Such pits form by the impact of an iceberg keel with the seabed Woodworth-Lynas et al., 1985), with the adjacent asymmetrical berms of sediment (Fig. 14E) likely having been excavated and deposited to one side by the iceberg upon impact. The connection of several of the depressions by curvilinear scours (Fig. 14B and D), marking the iceberg pathway across the seafloor in between low tides, strongly supports genesis by iceberg grounding (Brown et al., 2017;Newton et al., 2016). Furthermore, the depressions within the Anu survey were observed basinward of a major grounding zone wedge ( Fig. 14B and C), marking a potential ice sheet calving margin . Icebergs reworking the seafloor have been observed in water depths up to, and occasionally exceeding, 1 km offshore Greenland, e.g. Kuijpers et al. (2007), meaning that icebergs impacting the seafloor in the observed water depths of ~650 m is not uncommon.

Geothermal Gradient and Heat Flow
Geothermal gradient and heat flow estimations from this study provide the only data points across the entire continental shelf of western Greenland (Lucazeau, 2019). An average shallow geothermal gradient of 49 • C/km and heat flow of 70-90 mW/m 2 was estimated (Fig. 13), with these values being comparatively high when compared to the closest available (~225 km away from the study area) temperature probe data from within the centre of Baffin Bay (but only for the top <~2 m of sediment), which recorded five geothermal gradient values between 26 and 48 • C/km and heat flow values between 54 and 64 mW/ m 2 (Lucazeau, 2019;Pye and Hyndman, 1972). The heat flow estimate is also slightly greater than the nearest points (64-75 mW/m 2 ) in the global heat flow model by Davies (2013), albeit at the edge of its interpolated coverage.
Several factors may have resulted in higher geothermal gradient and heat flow estimations (Fig. 13). This could be associated with uncertainty in both the BSR based estimation method (up to ~50-60%; Grevemeyer and Villinger (2001); see section 3.3) and the temperature probe measurements of which the data are compared (Pye and Hyndman, 1972). An increase in heat flow is expected however, in areas of focussed fluid flow (Ganguly et al., 2000;Minshull and Keddie, 2010), such as above the MBR, and also due to the advection of heat towards underlying paleo-rift topography (the MBR) (Serié et al., 2017). Hydrocarbon leakage from the underlying Eocene gas reservoir may also add to this effect (Figs. 8-9 and 12). Heat flow may have been over estimated however, due to the use of seismic interval velocities to derive thermal conductivity (see section 3.3). The velocity within the GHSZ has likely been elevated by the presence of gas hydrate (Fig. 4) (Shankar and Riedel, 2011;Stoll et al., 1971), as such rapid lateral variations in velocity are unlikely to be geological. The thermal conductivity of methane hydrate is approximately equal to the pore water (<10% difference) (Ruppel, 2000;Waite et al., 2009), so elevated velocity would cause an over estimation of both thermal conductivity and heat flow, and the results may in fact be closer to those estimated previously (Lucazeau, 2019;Pye and Hyndman, 1972). Additionally, heat flow may not be in a Fig. 15. A composite seismic cross section in two-way-time across the Anu and Pitu 3D seismic surveys, overlain by an interpretation of the structure, stratigraphy and complex fluid migration history that likely characterises the study area. Fluid migration into the Cenozoic stratigraphy is observed via syn-rift carrier beds and tectonic faults from both within the ridge and above its western flank. The location of shallow gas accumulations, the Eocene gas reservoir and gas hydrate deposits, along with the potential migration pathways for these fluids is shown, along with possible effects of glaciation on hydrocarbon remigration. The location of the seismic line is shown on Fig. 3. steady state, due to the impact of recent glaciations on the temperature regime (Johansen et al., 1996), such as is observed on the Norwegian margin (Jung and Vogt, 2004;Mienert et al., 2000). Therefore, the hydrate system may not yet have fully re-equilibrated post-glaciation, affecting the BSR depth and heat flow estimates (see section 5.4).
Several localised, high geothermal gradient and heat flow anomalies are observed across the BSR area (Fig. 13), but no seismic features, that may suggest focussed fluid flow pathways (e.g. pipes or faults), are observed beneath (Ganguly et al., 2000;Minshull and Keddie, 2010). Instead, many of the anomalies coincide with either positive, likely depositional, or negative, erosional seabed features that likely formed during the last glaciation (Figs. 13-14) . These include multiple concentric and linear anomalies that represent iceberg pit-marks, glacial lineations and shallow seabed depressions (Fig. 13). These erosional features have caused localised deepening of the seabed, which is not reflected at the BSR depth (or at least is not visible given the seismic resolution), causing a thinner GHSZ, and increased apparent geothermal gradients ( Fig. 13A and E).

Possible Impacts of Glaciation
Melville Bay has experienced multiple cycles of ice sheet advance and retreat across the shelf since ~2.7 Ma , resulting in the removal of vast amounts of sediment from the inner shelf and the re-deposition of this sediment to the shelf edge, as part of the Melville Bay Trough Mouth Fan (mu-A) (Figs. 1-2, 11 and 15) Newton et al., 2020). These processes have transformed the structural position of the stratigraphy across the shelf (Fig. 2), with structural changes likely exacerbated by isostatic compensation resulting from load redistribution (Fjeldskaar and Amantov, 2018;Zieba and Grover, 2016). The combination of erosion from the inner shelf and the increased burial (and compaction) of post-rift sediments beneath the glacigenic wedge on the outer shelf (Fig. 2), as well as isostatic compensation, has led to the Melville Bay margin tilting basinward west of the MBR. Tilting of the stratigraphy has likely effected fluid migration by focussing fluids up-dip towards the inner shelf and above the MBR, possibly contributing to increased hydrocarbon occurrences in this area (Figs. 2, 6-7 and 15).
Long-term sediment redistribution and tilting likely occurred alongside repeated cycles of glacial loading on the shelf, impacting both isostasy and subsurface conditions (Cavanagh et al., 2006;Medvedev et al., 2019;Ostanin et al., 2017). Cyclic fluctuations in subsurface pressure and temperature during glacial-interglacial cycles, may have promoted fluid remigration (Fig. 15) (Goffey et al., 2016;Plaza-Faverola et al., 2011), with increased pore pressure beneath the ice load, potentially causing the fracturing of hydrocarbon seals and leakage (Hermanrud and Nordgard Bolas, 2002;Tasianas et al., 2018). This process may explain the dense network of fractures and gas pockets observed within the Eocene gas reservoir seal (mu-D1) and the overlying stratigraphy (Figs. 8-9). Furthermore, faults may reactivate, increasing fault permeability during deglaciation and load removal (Løtveit et al., 2011;Ostanin et al., 2017), possibly promoting leakage through the mu-D1 fractures, and increasing fluid migration through larger syn-and postrift faults parallel to the MBR. This may have contributed to the increase in hydrocarbon occurrences (Figs. 6-8 and 12) and the formation of gas hydrates (Figs. 10 and 15) above the western flank of the MBR.
Glacial-related variations in subsurface conditions can also affect the stability of gas hydrate (Portnov et al., 2016). The relatively large water depths in the study area (~650 m) mean the seabed is well within the stable region of the gas hydrate phase diagram (Fig. 16) and can thus be considered relatively stable, with hydrate dissociation at the seabed only likely to occur in response to glacial erosion (unlike that seen at shallower water depths in the Barents Sea, e.g. Andreassen et al., 2017).
These variations are more likely to have affected hydrate formation and dissociation at the phase boundary at the base of the GHSZ (Fig. 16) (Ruppel and Kessler, 2017;Serov et al., 2017). Ice loading and erosion during glacial periods likely cooled the underlying sediments and increased pore pressure, allowing increased hydrate stability and a deeper phase boundary than during ice-free periods (Fig. 16) (Grassmann et al., 2010;Ostanin et al., 2017). This suggests gas hydrate may have expanded downwards during glacial periods and once deglaciation and load removal began, would have experienced hydrate dissociation at the base of the stability zone (Portnov et al., 2016;Serov et al., 2017). This process may have contributed to the numerous free gas accumulations that are observed trapped at the base of the present day GHSZ (Figs. 6-10).
Post-glacial re-equilibration and shallowing of the phase boundary is evidenced beneath a glacigenic wedge (Figs. 10A and 13), where the BSR has shallowed significantly and mirrors the post-glacial seabed topography. However, evidence for the phase boundary not yet having re-equilibrated post-glaciation is observed below an area of recent glacial-related seabed erosion (evidenced by glacial lineations) ( Fig. 13C and 14) (Milkov and Sassen, 2000;Wang et al., 2006). This is causing an un-expectedly thin GHSZ in an area of deeper water ( Fig. 13C and E) and has likely increased geothermal gradient and heat flow estimates (Fig.  13).
The sensitivity of the gas hydrate deposits to future climate warming is an important consideration, especially within the Arctic, which is an area sensitive to the amplification of temperature increases compared to lower latitudes (Serreze et al., 2000;Serreze and Barry, 2011). Bottom water temperatures in Baffin Bay are estimated to increase by 0.2 to 3 • C by 2100 in a response to global warming (Stocker et al., 2014). However, due to the hydrate stability imparted by the relatively deep, and cold waters in Melville Bay, the gas hydrate deposits documented here are, at least in the short-term, likely less sensitive to projected future warming scenarios than hydrates in shallower waters, such as the Barents Sea . Even at the maximum of the estimated warming (3 • C warmer ocean bottom temperatures), hydrate at the seabed is likely to remain stable (Fig. 16). This scenario would Fig. 16. Temperature pressure conditions and methane hydrate stability. The phase boundary represents the transition of stable gas hydrate to free gas via dissociation and is based on the empirical relationship between pressure and dissociation temperature from Dickens and Quinby-Hunt (1994). Present day hydrate stability is shown at a seabed depth of 650 m and the present water temperature of 1.5 • C, with the estimated geothermal gradient through the gas hydrate stability zone (GHSZ) indicating the depth of the phase boundary and the observed bottom simulating reflector (BSR). Gas hydrate stability at the seabed is also shown to be less sensitive to future climate warming scenarios (3 • C increase in water bottom temperature in Baffin Bay) due to the relatively large water depths.
cause the hydrate phase boundary to shallow by ~50 m causing further dissociation at the base GHSZ. But for gas hydrate to dissociate at the seabed and be released into the atmosphere, ocean bottom water temperature would likely need to increase by ~8 • C, a scenario that is currently unlikely.

Conclusions
Numerous fluid related anomalies have been identified and mapped throughout the Cenozoic sedimentary succession of the Melville Bay continental shelf margin, using extensive 3D seismic data coverage. These data provide the most widespread and high-resolution imaging of the subsurface along the west coast of Greenland and the identified fluid anomalies represent the most northerly recorded occurrence of shallow gas and gas hydrates in Baffin Bay. Shallow gas anomalies have been categorized by different trapping styles, with free gas occurring in isolated pockets (category I), trapped at the base of the gas hydrate stability zone (GHSZ) (II), within laterally extensive glacial clinoforms (III) and finally within an Eocene aged, giant mass transport deposit reservoir (IV). Extensive gas hydrate deposits have been identified across an area of 537 km 2 via the identification of a discontinuous bottom simulating reflector (BSR) that represents the base of the GHSZ. The BSR was used to estimate a near surface geothermal gradient of 49 • C/km across the GHSZ and a heat flow of 70-90 mW/m 2 , providing the first geothermal gradient and heat flow data points on the entire west Greenland margin.
Paleo-rift topography is predicted to have influenced fluid migration and the thermal regime, with increased concentrations of gas and all gas hydrate deposits occurring above the centre and western flank of the Melville Bay Ridge (MBR). Regional post-rift stratigraphy onlaps and is tilted towards the ridge, promoting fluid migration towards the ridge crest. Related sets of faults along the ridge flanks have likely provided fluid migration pathways from potential deeply buried Cretaceous source rocks towards the overlying Cenozoic stratigraphy. Leakage from the Eocene reservoir is also likely, evidenced by condensed gas anomalies both within and above the densely faulted overlying seal.
Multiple cycles of glaciation since ~2.7 Ma have likely affected the contemporary distribution of hydrocarbons and the fluid migration history. Extensive sediment redistribution from the inner to outer shelf, along with the resulting isostatic compensation, has caused a regional tilt of the margin basinward, focussing fluids up-dip towards the inner shelf and above the ridge. Variations in subsurface conditions due to repeated cycles of glacial (un)loading, have likely promoted fluid flow, potentially causing seal breach and leakage from the Eocene reservoir due to increased pore pressures, as well as increased migration through ridge flanking faults due to fault reactivation. Finally, variations in pressure and temperature likely led to the expansion of gas hydrate deposits during glacial periods, and subsequent hydrate dissociation at the phase boundary during deglaciation; a process which may have contributed to the numerous free gas accumulations observed at the base GHSZ at present day. The relatively large water depths in this gas hydrate province would probably ensure continued stability of the top of the GHSZ both during de-glaciation in the past and during most future warming scenarios over the next several decades.

Data Availability
The seismic data used within this study were kindly provided by Cairn Energy PLC (Pitu), Shell (Anu) and TGS (Regional 2D). To access these data, the respective company or the Greenland authorities (MLSA) will need to be contacted directly.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.