The regional average geothermal gradient of 25.5°C/km is derived from the three BSR depths where sediment temperature is assumed to be at steady-state: away from the heat-conductive shallow salt, outside of the slide escarpment and outside of the MTD (red, yellow and white stars in the inset of Figure 2A). The 2D heat flow model shows a pronounced heat flow increase directly above the top of salt and elevated heat flow ~4000 m southward and northward along the profile c-d (Figure 3, see methods). The geothermal gradient predicted by the 2D heat flow model over the salt body is 30°C/km, and it gradually decreases to the average regional 25.5°C/km ~4000 m northward from the salt summit (Figure 3). However, these variations are insufficient to explain the shape of the irregular BSR, which requires geothermal gradients increasing from 16.3°C/km downslope to as high as 43°C/km above salt (Figure 2A).

Moreover, the 2D heat flow model predicts where the base of GHSZ would be at steady-state, and indicates that the sediments are still undergoing the residual post-slide temperature adjustment (Figure 3). To analyze this adjustment, profile c-d can be divided into three areas (Figure 2B): 1) the upslope area where the removal of the overburden is cooling the shallow sediments and drives the base of GHSZ down (Figure 2B, Figure 3), 2) the downslope area where the warming effect due to the deposition of the MTD drives the base of GHSZ upward, and 3) the areas outside of the MTD and slide escarpment where the GHSZ is relatively steady-state (Figure 2B, Figure 3). The fact that the BSR still closely mimics the pre-slide seafloor implies that the temperature adjustment is not yet significant, and that the landslide is relatively young. The following time-transient heat flow model provides constraints on the age of the slope failure.

**1D Transient Heat Flow Models**

We model the transient post-slide temperature field at the upslope slide location where sediments experience cooling, and at the downslope MTD location, where the sediments undergo warming (Figure 2A, B, see methods). We use a bottom water temperature of 4.2°C for the upper boundary condition. To constrain the lower boundary at each location, we apply constant heat flow that corresponds to the steady-state geothermal gradients predicted by the 2D heat flow model: 30°C/km at the upslope location and 25.5°C/km at the downslope location (Figure 3). Based on the 1D heat flow modeling (see equation (1) in methods), we define the transient temperature changes in the sediment column and derive the temperature profiles at certain times (time-temperature profiles) after the slide for both locations (see methods). Finally, to define the age of the slide event in the Orca Basin, we find the crossover of the three functions: the methane phase boundary curve, the observed BSR depth (reflecting the modern temperature), and the corresponding time-temperature profile (insets of Figure 2B).

2D steady-state conductive heat flow model used to estimate the effect of the salt body on the geothermal gradients along seismic section c-d (located in Figure 1). The model domain is 14 km deep and 16.5 km wide with constant basal heat flow. The properties of the sediment (i.e. thermal conductivity) vary with porosity and the properties of the salt vary with temperature (see methods and supplementary material). The model predicts an elevated geothermal gradient over salt (~30°C/km) and regional average geothermal gradient below the MTD (~25.5°C/km), which is however insufficient to explain the observed shift in the BSR in the Orca Basin (Figure 2A). We use the observed gradients at the upslope (30°C/km) and downslope (25.5°C/km) locations for our one-dimensional simulations (Figure 4).

**Upslope Location**

Given the 30°C/km steady-state geothermal gradient, the pre-slide temperature at the level of the modern seafloor (~220 meters below the pre-slide seafloor) was 10.7°C (Figure 4A). After the instantaneous ~220-meter seafloor drop caused by the slide, these warmer sediments were exposed to cooler bottom waters with a temperature 4.2°C (Figure 4A). Over thousands of years, the temperature within the subseafloor sediments gradually cools to adjust to the new boundary condition as shown by the time-temperature profiles (Figure 4A).

Based on the steady-state geothermal gradient, the pre-slide BSR would have been at ~ 465 m below the pre-slide seafloor (~245 meters relatively to the modern seafloor) (Figure 4a, black 0 kyr profile), and it will reach its post-slide steady-state depth at ~505 mbsf approximately 200 kyr after the slide event. Time-temperature profiles show the base of the GHSZ approaches approximately half-way to its steady-state depth during the first ~20 kyrs after the slope failure. Figure 4a shows that the intersection of the modern BSR (~342 mbsf, blue line in Figure 4a) and methane phase boundary curve corresponds to the ~8 kyr time-temperature profile (red curve in Figure 4A), which constrains the upslope location age estimate for the Orca submarine slide.

**Downslope Location**

At the downslope location, we model a 400-meter thick sediment mass added to the top of the pre-slide seafloor to simulate the deposition of the MTD (Figure 4B). The initial temperature profile within the MTD is somewhat ambiguous because slope failure is a chaotic process involving sediment redeposition and unpredictable rates of mixing with the cold bottom water. We select a uniform 7.5°C temperature throughout the MTD for the model, which was the mean temperature in the ~220 m thick upslope sediment column prior to failure (Figure 2B, 4B). The temperature at the top of the MTD is 4.2°C.

The model shows that at the downslope location the depth of the pre-slide BSR was ~720 meters below the pre-slide seafloor (~1120 mbsf) (Figure 4B). It will reach its post-slide steady-state depth at 680 mbsf approximately 250 kyrs after the slide event, whereas complete steady state will be reached in ~350 kyrs (Figure 4B). The modern BSR is observed in the seismic data at ~1090 mbsf, which is only 30 m above its pre-slide location (Figures 3, 4B). The modern BSR depth and methane phase boundary intersection corresponds to the ~14 kyr time-temperature profile (red line in Figure 4B), which constrains the downslope location age estimate for the Orca submarine slide.

A) Transient 1D heat flow modeling at the upslope location (Figure 2B) showing the time-temperature profiles after the slope failure (removal of overburden). The figure shows the intersection of methane phase boundary (black curve) and modern BSR depth (blue line) corresponds to an 8 kyr time-temperature profile (red line) indicating the age of the Orca slide derived from the upslope location. B) Transient 1D heat flow modeling at the downslope location (Figure 2B) showing the time-temperature profiles after the MTD deposition. The model shows that the intersection of methane phase boundary (black curve) and modern BSR depth (blue line) corresponds to a 14 kyr time-temperature profile (red line). For the upper boundary condition, we use a bottom water temperature of 4.2°C. At the lower boundary, we apply constant heat flow that corresponds to 30°C/km geothermal gradient at the upslope location and 25.5°C/km at the downslope location (see supplementary material for thermal conductivity and diffusivity depth profiles). Front panel shows that both slide ages determined for the Orca Basin correlate with a rapid post-glacial sea level rise40, which could cause increased slope instability across the worldwide continental margins13.

**Age Of The Orca Submarine Slide**

At the upslope location, we predict an age of ~8 ka for the Orca submarine slide, whereas at the downslope location we predict an age of ~14 ka. We consider the 8 ka age estimate to be more accurate for two reasons. First, sediment removal upslope was likely a single fast-moving event. In contrast, the MTD at the downslope location may be an amalgamation of several landslides that occurred around the rim of the Orca basin 33. Thus, the temperature profile within the sediments at the downslope may record a number of slides, some older than the one released at the upslope location. Second, the depth interval between the pre-slide and modern BSR is wider and shallower at the upslope location compared to the downslope location (Figure 4A, B). Within this interval upslope, the thermal signal propagates faster resulting in the wider-spaced time-temperature profiles. Thus, the upslope location provides better age resolution than the densely clustered profiles at the downslope location.

The impact of slope failures on gas hydrate accumulations is different in the interval below the MTD than below the slide escarpments. Below the MTDs, there would be gas hydrate dissociation as the base of GHSZ rises with the temperature increase (Figure 2B). Gas hydrate dissociation below the MTDs (i.e. the downslope location) is an endothermic reaction41 that would draw the heat from the sediments, slow the general warming trend, and consequently increase the modeled slide age. In contrast, below the slide escarpments (i.e. the upslope location), the deepening GHSZ entraps the underlying gas, which forms hydrate (Figure 2B). Gas hydrate formation is an exothermic process accompanied by heat release, which may slow the general cooling trend and likewise increase the modeled slide age. Therefore, the ~8 ka is the youngest age estimate for the Orca landslide. Accurate quantitative characterization of hydrate formation and dissociation depends on gas hydrate and gas saturations over the area, which are not available. In any case, our modeling shows that the effect of slope failures on gas hydrate systems can continue over tens of thousands of years after the slide events.

Our inferred age of the Orca slide (~8 ka) coincides with the final stage of rapid late Pleistocene–early Holocene sea level rise in the Gulf of Mexico averaging at ~12 mm/yr between 14.5 and 8 ka40(inset of Figure 4B). An earlier study showed that this fast post-glacial sea-level rise could induce seismicity along faults on passive continental margins and increase the frequency of submarine slides worldwide13 (inset of Figure 4). A similar mechanism could be responsible for the activation of deep-rooted and supra-salt faults, triggering submarine slides in the Gulf of Mexico, including the Orca Basin.

**Model Sensitivity**

Our landslide age estimates may be affected by factors that control gas hydrate stability, such as the presence of heavy hydrocarbons (ethane, propane, butane, pentane), and high salinity. In addition, geologic factors can also influence hydrate stability by altering temperature and pressure, such as elevated post-slide sedimentation rates and eustatic sea level change. These factors can contribute to the depth of the GHSZ, yet the impact on our age estimates is unlikely to be significant.

First, the presence of heavy hydrocarbons and/or elevated salinity in the shallow subseafloor sediments would provide the opposite effects at the upslope and downslope locations. Higher salinity would shift the methane hydrate phase boundary leftward in Figure 4 and result in a younger modeled slide age at the downslope location and older slide age at the upslope location. The presence of heavy hydrocarbons would shift the methane-hydrate phase boundary rightward in Figure 4 and result in the older modeled slide age at the downslope locations and younger age upslope. Neither salinity anomalies nor thermogenic gas are supported by our well log and seismic data.

Second, the high post-slide sediment deposition rate would result in a younger modeled slide age at the downslope location and older age at the upslope location. Based on the 3D seismic data, the hemipelagic drape within the slide escarpment is absent or <8 m thick, indicating that the effect of the sediment deposition is negligible. Additionally, the thin layer of post-slide sediments supports the relatively young age of the slide.

Third, the post-glacial sea level rise in the Gulf of Mexico40 would only result in <7 m uncertainty in the pre-slide base of GHSZ based on the corresponding shift of the methane-hydrate phase boundary profile. This slightly affects the pre-slide BSR location but provides negligible effect on the slide age.

**Universal Approach To Submarine Slide Dating Using BSRs**

Equation (2) (see methods) provides a simple analytical method for a quick-look slide age (\({t}_{slide}\), s) estimate using the modern BSR depth (\({z}_{bsr}\), m) below slide escarpments (similar to the Orca upslope location) with a known temperature at the BSR depth (\({T}_{bsr}\))(Figure 2A). We further find that parameter \(\frac{{z}_{bsr}-{z}_{bsr,o}}{{z}_{bsr,1}-{z}_{bsr,o}}\) can be used to quantify the fractional heat dissipation after the slide, with a value of zero referring to the initial condition immediately after the slide and a value of 1 referring to the post-slide steady state (Figure 4A, 5). The plots of the fractional heat dissipation (\(\frac{{z}_{bsr}-{z}_{bsr,o}}{{z}_{bsr,1}-{z}_{bsr,o}}\) ) versus the dimensionless time (\(\frac{{\kappa }t}{{z}^{2}bsr}\)) all fall into one curve at different locations with different water depths, temperature gradients, thermal diffusivities, BSR temperatures and/or landslide thicknesses (Figure 5: red curve). This means that the age of a submarine landslide \({t}_{slide}\)can be estimated simply from the diagram shown in Figure 5 with easily accessible parameters (\({z}_{bsr}\), \({z}_{bsr,o}\), \({z}_{bsr,1},\) \({\kappa }\)). It is evident that the BSR shift from the pre-slide to post-slide steady-state \(({z}_{bsr,1}-{z}_{bsr,o}\)) can be roughly equated to the thickness of a submarine slide (Figure 4A, 5).

The analytical solution (see equation (2) in methods), which estimates submarine slide age using BSR depths under the landslide escarpments, similar to the upslope location at Orca. \({z}_{bsr}\) is calculated from equation (2) by assuming a constant, known temperature at the BSR depth (\({T}_{bsr}\)) during the evolution. \({z}_{bsr,0}\) is the depth of a pre-slide steady-state BSR that can be estimated based on the thickness of removed sediments; \({z}_{bsr,1}\) is the depth of a post-slide steady-state BSR; κ is the average thermal diffusivity of the subseafloor sediment (see supplementary material); \({t}_{slide}\) is the age of a submarine landslide. The analytical solution for the Orca upslope location using 30°C/km geothermal gradient (red dashed arrows) dates the slide to ~7.8 ka, which is similar to the result of our numerical simulation. Gray dashed lines show the sensitivity of the analytical solution to the varying geothermal gradients.

The age of the Orca slide from the analytical solution is ~7.8 ka, which is in a good agreement with our numerical solution (Figure 5). The discrepancy between the analytical and numerical solutions may be caused by a difference in domains; the numerical solution is solved for a finite domain, whereas the analytical solution is solved for a semi-infinite domain. There are also slightly different physical properties with depth that could cause variation in the results as the numerical solution has varying porosity and thermal diffusivity, whereas the analytical solution assumes homogeneous sediments properties with depth.

Our seismic and modeling-based landslide dating technique is a novel method for dating submarine landslides without core data. Because gas hydrate systems occur worldwide and often coexist with submarine landslides, there may be many locations where our approach can be applied. Moreover, this technique is especially relevant given the expanding public seismic databases worldwide. The following are only few examples where published seismic data appear to show irregular BSRs deviating from the seafloor bathymetry below the landslide escarpment and/or below MTDs: the Cape Fear slide complex offshore the US East Coast 42, offshore Oregon, USA 43, the Storegga slide offshore southern Norway44, the Brunei slide offshore Brunei45, and the Hinlopen megaslide offshore Svalbard46,47.