Establishing how past climate change affected the stability of marine methane hydrate is important for our understanding of the impact of a future warmer world. As oceans shallow toward continental margins, the base of the hydrate stability zone also shallows, and this delineates the feather edge of marine methane hydrate. It is in these rarely documented settings that the base of the hydrate stability zone intersects the seabed and hydrate can crop out where it is close to being unstable and most susceptible to dissociation due to ocean warming. We show evidence for a seismically defined outcrop zone intersecting canyons on a canyon-incised margin offshore of Mauritania. We propose that climatic, and hence ocean, warming since the Last Glacial Maximum as well as lateral canyon migration, cutting, and filling caused multiple shifts of the hydrate stability field, and therefore hydrate instability and likely methane release into the ocean. This is particularly significant because the outcrop zone is longer on canyon-incised margins than on less bathymetrically complex submarine slopes. We propose considerably more hydrate will dissociate in these settings during future ocean warming, releasing methane into the world’s oceans.

Naturally occurring methane hydrate is abundant in oceanic settings, but changes in pressure and temperature (P-T) cause it to dissociate, potentially allowing methane to escape to the atmosphere and amplify climatic change (Ruppel and Kessler, 2017) or oxidize in seawater and contribute to ocean acidification (Biastoch et al., 2011). Seaward shifts of the updip feather edge of marine methane hydrate resulting in methane venting to the ocean have been documented due to contemporary warming of bottom waters (Skarke et al., 2014; Ketzer et al., 2020), but this has not been demonstrated for past climate change on canyon-incised margins. There are >5800 submarine canyons, with a total length of 254,000 km, on Earth (Harris and Whiteway, 2011), and methane emissions and hydrate have been detected in canyons (e.g., Barkley Canyon, offshore British Columbia, Canada: Thomsen et al., 2012; Hudson Canyon, offshore northeastern United States: Weinstein et al., 2016). But a bottom-simulating reflection (BSR) marking the base of the hydrate stability zone (BHSZ) has not been identified intersecting canyon floor and wall settings where the vulnerable hydrate could crop out. Determining whether there has been any movement of the feather edge due to past ocean warming, the evidence for this, and the relative importance of canyon-incised slopes as opposed to less complex settings requires consideration so that we can better understand the role this setting will play during future climatic change.

Petroleum industry three-dimensional (3-D) seismic data have imaged the head of the Cap Timiris Canyon system and contourite moats (Fig. 1) on the Mauritanian margin. The feather edge of the marine methane hydrate and outcrop zone have been described to the southeast (Davies et al., 2015) (Fig. 1C). The petroleum industry borehole Ras Al Baida-A1 (location in Fig. 1A) shows that the succession is fine-grained, Quaternary in age, and during this time the sedimentation rate was ~19.2 cm k.y. −1 (Dahi et al., 2013).

The 3-D data (Figs. 1, 2, and 3) were acquired in 2012 and have a vertical resolution of ~8 m. The seismic volume was depth converted and is zero-phased. An increase in acoustic impedance is shown as a red-black-red reflection. The BSR is a reflection with the opposite polarity to the seabed or a series of aligned reflection terminations (Fig. 1F) (c.f. Minshull et al., 2020). We also identify relict BSRs, which are probably due to residual gas at the location of an earlier position of the BHSZ (Zander et al., 2017) on the basis of its cross-cutting relationship with stratal reflections and parallel relationship to the modern BSR. Below the BSR, hydrate is not stable and free gas occurs (the free-gas zone), and it is this juxtaposition that causes the reflection. We used two methods to determine whether there have been shifts in the BHSZ. Firstly, we identified lines of intersection (LoIs) using root mean square (RMS) amplitude maps of seismic reflections. The LoIs mark the intersections of the reflections with the earlier and modern positions of the BHSZ manifested by changes in the seismic amplitude of the stratal reflection. These changes in amplitude are linear and mainly parallel to each other in plan view. LoIs can occur within the HSZ or free-gas zone and can indicate multiple resets of the BHSZ that can be accounted for by past changes in the P-T conditions of bottom waters (Davies et al., 2017). Secondly, the relict BSRs show the BHSZ has reset. We modeled the past positions of the BHSZ to compare with the relict BSRs. The modeling method is exactly the same as that of Davies et al. (2017), but here we used eight climate models rather than one, to test for a match with the seismic observations of the LoIs and relict BSRs (CCSM4, GISS-E2-R, MPI-ESM-P, IPSL-CM54-LR, CNRM-CM5, MIROC-ESM, MRI-CGCM3, and FGOALS; URLs for the models are provided in Table S1 in the Supplemental Material1). Changes in relative sea level were from a time series by Bintanja and van de Wal (2008). The modern BHSZ was calculated using a vertical profile of measured modern ocean temperatures characteristic of the wider oceanic area around the Cap Timiris Canyon. The profile was extracted from the World Ocean Atlas (https://www.nodc.noaa.gov/OC5/indprod.html; Locarnini et al., 2013) and consisted of annual mean temperatures averaged over the region 18°–22°N, 21°–17°W. Changes in the temperature profile of the sediment due to variations in bottom-water temperature were calculated vertically using a one-dimensional uniform and constant heat diffusivity of 10−6 m2 s−1 (e.g., Muraoka et al., 2014) and a geothermal gradient that is free to evolve within the sediments in response to temperature changes at the bottom of the water column, but is prescribed to be equal to 26 °C km−1 10 km below the seafloor (Macgregor, 2020). Given typical hydrate volumetric concentrations of 1%–10%, heat fluxes associated with the formation or dissociation of hydrates would be 10×–100× smaller than typical marine geothermal heat flows of 0.1 W m−2 (e.g., Hofmann and Morales Maqueda, 2009) and have therefore been ignored.

There are 16 canyons in the Cap Timiris Canyon system, which are 0.05–0.4 km deep and 0.1–3 km wide (Figs. 1A1F). For each canyon, the northwestern canyon wall differs markedly from the southeastern one. A typical northwestern canyon wall (Fig. 1E) comprises a series of stacked, buried seismic packages (e.g., marked 1–5, Fig. 1E). In Figure 1E, each is ~0.5–1 km wide, and they step progressively from northwest to southeast. In contrast, reflections on the opposite canyon wall show none of the packages, and instead, reflections are truncated (marked “T” in Fig. 1E).

The BSR progressively shallows in a landward direction and intersects the seabed at a water depth of ~720–730 m (the seabed intersection depth) (Fig. 1F). The 720–730 m isobath is deflected landward by as much as 4.5 km due to the bathymetric complexity of the canyon-incised margin (green line in Fig. 1B), and where the canyons are best developed, the 720–730 m isobath is 199 km long. Without the lateral deflections caused by the canyons, the isobath would be ~84 km long (Fig. 1B, green line). Canyons can develop fractal and dendritic patterns, therefore 199 km is probably a conservative measurement for the length of the seabed intersection depth, where hydrate nearby is expected to be to be unstable.

In seismic lines that cross the canyons, the BSR is manifested as a series of aligned amplitude terminations that deepen, tracking the canyon walls (Figs. 2A and 2C). We mapped a number of stratal reflections between the canyons that intersect the downward-deflecting BSR and selected two, reflections X and Y, to test whether there are linear changes in seismic amplitude, consistent with the occurrence of a LoI that could be interpreted as evidence for shifts through time in the position of the BHSZ (Figs. 2A and 2C). Maps of the RMS amplitude of both reflections (Figs. 2B and 2D) reveal a series of curvilinear changes in seismic amplitude with consistent spacing of 250–500 m (marked on Figs. 2B and 2D). In both examples, the features are on both sides of the LoI that marks the modern intersection of the BSR with each reflection (marked 1–5 in Figs. 2B and 2D) and they are parallel to it and the canyon wall. Seismic cross sections (Figs. 2A and 2C) confirm that the LoIs correspond to subtle increases in seismic amplitude (marked 1–5 in Figs. 2A and 2C). There is no evidence for relict BSRs in these examples.

At another canyon, where the seismic line follows the canyon margin rather than being orthogonal to it, the angle of intersection between the stratal reflection and the modern BSR is lower. The RMS amplitude map for reflection Z shows that there are two LoIs, one marking the modern intersection of the BSR and the reflection, and a shallower one located within the HSZ (Figs. 3A and 3B). In this case, there is also evidence for a throughgoing, deeper relict BSR (Fig. 3A). In Figure 3A, the modern BSR intersects the seabed at 720–730 m water depth. Additional examples of LoIs within the study area are provided in Figure S2 in the Supplemental Material1.

The northwest-southeast–oriented contourite moats (Fig. 1B) are due to southeastward-flowing deep-water currents (Peña-Izquierdo et al., 2012). The stacked seismic packages that are found on the northwestern side of the canyon margins are similar to laterally inclined packages in the South China Sea (Zhu et al., 2010), interpreted to be clinoforms that formed as bottom currents swept muddy sediments across the slope, over the canyon rims, and into the canyons (c.f. Shanmugam, 2008). This causes the canyons to migrate laterally in one direction and explains the seismic packages and erosion on the opposite canyon wall. Sedimentation rates dropped during interglacial times (Henrich et al., 2010), leading to the dominance of deep-water erosion, therefore erosion of the southeastern canyon margins should have been more significant during interglacial than glacial periods.

The intersection of the BSR on the canyon walls and floor (Figs. 1F and 3A) at 720–730 m water depth indicates that hydrate is located at or immediately below the seabed (a “seismically defined” outcrop zone; Davies et al., 2015). The canyons also cause the BHSZ to deflect downward because the water in the canyon has a localized cooling effect on sediment surrounding and below the canyon (Portnov et al., 2016). This sets up cross-cutting relationships between the BSR and reflections X and Y (Figs. 2A and 2C). The linear changes in seismic amplitude of reflections X and Y (Figs. 2B and 2D) are very similar to those identified in other 3-D seismic data sets where reflections cut across a BSR, and like in some of these examples, the changes occur on either side of the modern BSR in both the HSZ and free-gas zone (Davies et al., 2017). Following previous interpretations, we propose the linear increases in seismic amplitude are the result of relict gas accumulations in low-permeability sediment that were sealed vertically by hydrate-clogged strata, where the gas has not yet dispersed (Zander et al., 2017). Changes in gas saturation could also be preserved in the form of changes in hydrate saturation within the HSZ (Davies et al., 2017), and this could explain the linear changes in seismic amplitude within the HSZ (marked 1 in Figs. 3A and 3B). In the free-gas zone and the HSZ, they are manifested as high-amplitude bands ~0.2–1.0 km wide. The downdip boundaries to these bands, where high amplitude changes to lower amplitude, could represent the bases of relict gas accumulations (Figs. 2B and 2D). Where the seismic line is approximately parallel to the canyon, the angular relationship between the BSR and the mapped reflection is lower, and therefore the LoIs are more widely spaced (Fig. 3B).

Like other examples of LoIs, we propose the LoIs formed due to repeated shifts in the BHSZ due to changes in P-T conditions in the past (Davies et al., 2017). We can rule out a sedimentary origin for the amplitude changes because the features occur only on either side of the modern BSR and not in any of the other reflections adjacent to the canyons (Figs. 2A and 2C).

To test the hypothesis that ocean warming since the Last Glacial Maximum (LGM) caused the shift from the relict to the modern BSR and/or the LoIs, we modeled the position of the BHSZ during the LGM for the seismic lines (Figs. 2A, 2C, and 3A). The sea level relative to present that we adopt for the LGM is −123.40 m (Bintanja and van de Wal, 2008). All but one of the climate models resulted in a deeper BHSZ than the modern BSR during the LGM; for the seismic line in Figure 3A, two models (CCSM4 and GISS-E2-R) showed a very close match, and CCSM4 is marked on the seismic line and lies on the relict BSR (Fig. 3A). The bottom temperatures predicted by CCSM4 for the LGM are between 2.2 °C and 2.8 °C cooler than at present. Therefore, the location of the relict BSR is consistent with the probable deeper position of the BHSZ during the LGM ~20 k.y. ago, and bottom-water warming has caused it to shallow by >~200 m to its present location. It cannot be ruled out that this match is coincidental and the deeper relict BSR represents an earlier Quaternary glaciation. But gas should disperse by diffusion over time, and therefore one would expect the feature to be the result of the most recent cooler period rather than an older one. The modeled position is significantly deeper than any of the LoIs, which are much closer to the modern BSR (Figs. 2A and 2C). The regular spacing of the LoIs and their situation parallel to the erosional canyon margin (Figs. 2B, 2D, and 3B) indicate a relationship between erosion and the formation of the LoIs. Furthermore, the occurrence of a LoI within the present HSZ indicates that most recently there has been an episode of cooling, which reset the BHSZ from this location to its present deeper and more landward position. We propose that erosion of the southeastern canyon margins and therefore cooling were caused by contourite currents, consistent with the Holocene interglacial (the past 11.7 k.y.) when sedimentation rates dropped (Henrich et al., 2010) and deep-water erosion processes became more dominant.

We have identified two causes of shifts of the hydrate stability field, and therefore hydrate instability and methane release into the ocean: (1) ocean warming drives significant shallowing of the BHSZ and outcrop zone, which would have destabilized hydrate at and immediately below the seabed, probably allowing methane to enter the water column; then (2) a less-significant subsequent deepening of the HSZ due to canyon margin erosion. Therefore, it is likely that migration, cutting, and filling of canyons can drive shifts in the position of the BHSZ and the flux of methane into and out of hydrate reservoirs; this is consistent with observations made elsewhere (Jin et al., 2020).

Our work is one of a small number of documented examples of movement of the feather edge of marine methane hydrate due to ocean warming in the past, most likely since the LGM. We propose that the volume of hydrate that would have been unstable is significantly higher on canyon-incised margins than less bathymetrically complex submarine slopes, and therefore these margins should have a proportionately more important role to play in releasing methane into oceans in a future warming world. Seismically imaged relict BSRs and LoIs allied to modeling the past positions of the BHSZ using climatic models is starting to show that we can potentially untangle how climate, sedimentation, and erosion, along with other factors, control the dynamics of these vast methane reservoirs.

We are grateful to Schlumberger Limited for access to Petrel software for the seismic interpretation, BP for access to the seismic data, while they licensed the study area, and six anonymous reviewers for their comments. The National Natural Science Foundation of China (grant 42006066) is acknowledged for funding Ang Li’s research.

1Supplemental Material. Table S1 (website URLs for the climate models), Figure S2 (additional examples of lines of intersection), and Figure S3 (uninterpreted versions of Figures 2 and 3). Please visit https://doi.org/10.1130/GEOL.S.14390660 to access the supplemental material, and contact editing@geosociety.org with any questions.
Gold Open Access: This paper is published under the terms of the CC-BY license.