Where and Why Do Submarine Canyons Remain Connected to the Shore During Sea‐Level Rise? Insights From Global Topographic Analysis and Bayesian Regression

The efficiency of sediment routing from land to the ocean depends on the position of submarine canyon heads with regard to terrestrial sediment sources. We aim to identify the main controls on whether a submarine canyon head remains connected to terrestrial sediment input during Holocene sea‐level rise. Globally, we identified 798 canyon heads that are currently located at the 120m‐depth contour (the Last Glacial Maximum shoreline) and 183 canyon heads that are connected to the shore (within a distance of 6 km) during the present‐day highstand. Regional hotspots of shore‐connected canyons are the Mediterranean active margin and the Pacific coast of Central and South America. We used 34 terrestrial and marine predictor variables to predict shore‐connected canyon occurrence using Bayesian regression. Our analysis shows that steep and narrow shelves facilitate canyon‐head connectivity to the shore. Moreover, shore‐connected canyons occur preferentially along active margins characterized by resistant bedrock and high river‐water discharge.

Based on a global compilation of submarine canyons, Harris and Whiteway (2011) showed that shelf-incising canyons prevail along the western, tectonically active margins of the Americas that are characterized by high sediment supply. Smith et al. (2017) focused on the West Coast of the United States and found no correlation between canyon occurrence and shelf gradient. Instead, their analysis showed that coarse sediment from onshore catchments exposing durable bedrock governs offshore canyon-head incision, modulated by wave focusing of canyon bathymetry (Smith et al., 2017(Smith et al., , 2018). Yet, whether onshore processes and lithological composition control headward submarine canyon erosion at the global scale has not been investigated.
We study global patterns of shore-connected canyons to identify the main controls on their occurrence. Our analysis is driven by two hypotheses: First, we hypothesize that submarine canyon heads remain connected to the shore upon postglacial sea-level rise if the shelf is narrow and steep. Second, we test the hypothesis of Smith et al. (2017Smith et al. ( , 2018) that submarine canyon heads remain connected to the shoreline when located offshore tectonically uplifting regions with durable bedrock. We categorize canyon heads as "shore-connected" and as "close to the 120m-depth contour" and assert that the latter may have been connected to the shore during the LGM but were unable to erode headwards to the present-day shoreline. We then predict shore-connected canyons using Bayesian penalized regression, a technique that selects predictor variables that are relevant for predicting a specific response variable.

Methods and Data
Our study relies on the Shuttle-Radar-Topography Mission (SRTM30_PLUS) 30-arc-second database (∼1 km resolution at the equator) (Becker et al., 2009). We excluded islands and oceanic plateaus and limited our analysis to 50°N and 50°S where predictor variables (Table 1) are available.

Submarine Canyon Variables
We used the data set from Harris et al. (2014) who mapped "large" canyons that extend over a depth range of at least 1,000 m and are incised at least 100 m. Hence, this analysis focuses on large features and omits smaller canyons. We manually mapped canyon heads, assigned them to the submarine-canyon polygon of Harris et al. (2014) (Figure 1), and computed the shortest Euclidean distance from each canyon head to the present-day shoreline and to the 120m-depth contour (LGM shoreline). Canyon heads located <6 km away from either shoreline were classified as "shore-connected present" or as "120m-contour canyons", respectively (Figures 1a and 1b). The distance exceeds the offshore limit of longshore sediment transport (up to 5 km; Sweet & Blum, 2016) by 1 km to account for uncertainties in mapping due to the latitudinally variable spatial resolution of the SRTM30_PLUS data set. We assume that "shore-connected canyons" receive terrigenous sediment directly from the shore through either longshore sediment transport or by direct river input. With this approach, we potentially include canyons that are not strictly connected to the ocean littoral cell but should capture all truly connected "large" canyons. We manually corrected canyon-head misclassifications during visual inspection. Although some canyons (Swatch-of-No-Ground, located 150 km off coast, Rogers et al., 2015;Indus canyon, 17 km, Li et al., 2018) presently receive terrestrial sediment through clinoform progradation, we did not classify these as shore-connected.

Terrestrial Variables
Topographic analysis was conducted with ArcGIS, MATLAB, and TopoToolbox (Schwanghart and Scherler, 2014). Terrestrial drainage-basin metrics (elevation, gradient, area, river steepness (ksn)), and river outlets were determined using the 500-m-grid and flow directions of the HydroSHEDS compilation (Lehner et al., 2008) (Table 1) The ratio of the 90th to 50th percentile of the mean annual rainfall dimensionless This ratio is a measure for how extreme rainfall is, high values indicate extreme rainfall events, area-weighted mean of the 90th/50th percentile of the mean annual rainfall for each catchment was assigned to each river outlet using the upslopestats function of TopoToolbox Boers et al. (2014) TRMM_90_50_weighted Extreme rainfall events may result in high instantaneous sediment supply and possibly in hyperpycnal flows flowing into canyon heads

Marine predictors All marine predictors are not weighted
Gradient of the adjacent shelf degree Computed outline of shelf shape file (Harris et al., 2014) We obtained estimates of water discharge (Qw) at each river outlet by integrating annual runoff from the Global Runoff Data Centre (Fekete et al., 2002). We used the BQART model for pre-human riverine suspended sediment flux (Qs) (Syvitski & Kettner, 2011), which plays an important role in offshore sediment transport (Hage et al., 2019). Bedload was estimated by two empirical equations (Table 1). We extracted a mean, area-weighted erodibility index for each catchment using the global erodibility index (GEroID) of Moosdorf et al. (2018), ranging from low-erodibility metamorphic rocks (=0.8) to high erodibility (=3.2) for unconsolidated sediments. Area-weighted means of annual rainfall and its variability were calculated using Tropical-Rainfall-Measurement-Mission data (TRMM, courtesy of Boers et al. [2014]). Finally, we assigned peak-ground accelerations (PGAs) to river outlets using the Global Seismic Hazard Map (Shedlock et al., 2003).

Weighting of Terrestrial Variables to Submarine Canyon Heads
To assign terrestrial variables to submarine canyons, we first determined the closest point on the shoreline (XY coast ) for each canyon head (Figure 1c), and then computed the distance between XY coast and river outlets on the adjacent continent. These distances i d together with the catchment areas i A of each river outlet i served as weights in a distance-weighted averaging approach. Specifically, we calculated the weights  i as: and assigned the weighted average to the corresponding canyon head. Our weighting scheme reflects that a one-to-one assignment of river outlets to canyon heads is often ambiguous and accounts for the increased importance of nearby and large catchment outlets ( Figure S3). Moreover, our approach exonerates us from choosing an arbitrary number of drainage basins that may shed sediment into a particular submarine BERNHARDT AND SCHWANGHART 10.1029/2020GL092234 6 of 15  canyon. We weighted the variables toward XY coast and not to the canyon head itself to avoid incorporating the canyon-head-to-shore distance into the predictor variables, as this is implicitly what we are aiming to predict.

Submarine Groundwater Discharge
Luijendijk et al. (2020) globally simulated coastal fresh submarine groundwater discharge (SGD). To assign weighted SGDs to canyon heads, we computed the centroid of each coastal watershed, assigned the SGD to that centroid, and weighted the value by its distance (d) to each XY coast of each canyon head using the weighting factor λ:

Marine Variables
To acquire the mean gradient of the continental slope in the vicinity of canyon heads, we extracted the outlines of the slope shapefile of Harris et al. (2014). We set elevation values within the extent of canyons to NoData and used a Laplacian interpolation to smoothly interpolate inward from these outlines. The technique is referred to as image inpainting (Stolle et al., 2019) and reconstructs a continental slope devoid of canyons. Analogously, we calculated the shelf gradient adjacent to each canyon head using the shelf shapefile of Harris et al. (2014). To calculate shelf width at each canyon head, we extracted pixels at the oceanward shelf boundary and calculated the shortest Euclidian distance to the shoreline. We chose the 100 nearest-neighbor pixels (∼86 km) along the oceanward boundary and calculated the mean shelf width from these pixels. We used a large number of nearest-neighbors to minimize the impact of canyon-head indentation into the shelf (Figure 1c). Results were inspected visually and corrected where wide canyons decreased shelf width. For the LGM shelf width, the same calculation was based on a digital elevation model (DEM) where 120 m of elevation where added to each pixel to simulate the LGM landscape. Shelf-edge depth was calculated by extracting the water depth from the 20 nearest-neighbor pixels (∼17 km) of the outer shelf boundary. We assigned storm-surge heights of 1-in-100 years extreme sea levels of Muis et al. (2016) to each canyon head. Bergsma and Almar (2020) extracted global wave heights and periods from ESA's Sentinel2 constellation and calculated the depth of closure (maximum water depth of littoral sediment transport; Hallermeier, 1980). We extracted wave height, period, and depths of closure for each XY coast and assigned these to adjacent canyon heads.

Predictive Modeling-Bayesian Penalized Regression
This study aims to identify the controls on continued shelf incision and maintenance of canyon-headto-shore connectivity. We computed a hexagonal net (50,000 km 2 /hexagon) and assigned the number of shore-connected canyon heads in each hexagon (Figure 1b). This is the number we predict (the "response"), using 34 predictor variables ("predictors," Table 1). To extract the weighted predictors for each hexagon, we computed the hexagon midpoints and their corresponding nearest location at the coast (XY coast _ hexgrid ) and weighted the predictors for individual canyon heads onto this coastal location (Equations 1 and 2).
To identify the most important predictors and to globally predict shore-connected canyon heads, we employed Bayesian penalized regression. Bayesian statistics apply probabilities to statistical problems offering a way to learn from new data to update prior beliefs while accounting for uncertainties (Efron, 2013;Korup, 2021). A frequentist approach to penalization is Lasso regression which uses a penalty term to shrink small regression coefficients to zero (hence reducing or eliminating unimportant predictors) (Tibshirani, 2011). In Bayesian penalized regression, penalization is incorporated through the choice of prior distribution (e.g., van Erp et al., 2019). We used bayesreg, a MATLAB toolbox for fitting Bayesian penalized regression models (Makalic & Schmidt, 2016). All predictors were centered and scaled. As we predict counts of shore-connected canyons per hexagon, we chose a Poisson distribution for the response (see Supporting Information [SI]). Our choice of shrinkage prior followed the procedure of van Erp (2020) (SI includes prior-sensitivity analyses, Figures S4-S7). Finally, we quantified the importance of each predictor adopting the Bayesian feature-ranking algorithm of Makalic and Schmidt (2011). The rank corresponds to the strength of the association between the predictor and the response (Tables S1-S5) and is based on the 75th percentile of the complete set of rankings for each posterior sample (SI includes a complete list of model parameters).

Shore-Connected Canyon Occurrence
Our data set comprises 4,633 canyon heads, of which 2,765 are classified as blind canyons and 1,702 as shelf-incising. From the latter, 798 were classified as 120 m-contour canyons and 183 as shore-connected canyons. 120m-contour canyons occur globally along passive and active margins (Figure 1a). In contrast, during today's sea-level highstand, most shore-connected canyons straddle along active margins (n = 114, Figure 2a) with spatial hotspots along the Mediterranean active margin, and the Pacific coast of central South America and Central America (Figure 1b). Moreover, shore-connected canyons occur frequently along the Californian coast, the Indian-Ocean coast of the Arabic Peninsula and the Eastern Black Sea. Isolated shore-connected canyons occur along the coasts of Africa (Figure 1b). Figure 2b shows the number of shore-connected canyons per hexagon plotted against the four highest-rank predictors (Figure 3d). 120m-contour canyons occur at shelf widths from <2 to 400 km ( Figure S2), whereas shore-connected canyons occupy narrow shelf widths from <2 to 31 km (Figure 2b). Only 8 shore-connected canyons occur at shelf widths between 20 and 31 km, 33 at shelf widths between 10 and 20 km, and the majority occurs at shelves <10 km wide. Shore-connected canyons occur preferentially where shelf gradients exceed 0.5°, the difference between the present-day and LGM shelf width is minimal (<27 km), and become most abundant at differences <9 km. One exception is the Congo canyon which occurs at a ∼50-km-wide shelf ( Figure S2). Shore-connected canyons occur along a wide range of erodibility indices but are absent at highly erodible catchment lithologies (GEroID > 2). The maximum number of shore-connected canyons (n = 20) is located offshore southern France.

Prediction of Shore-Connected Canyons
Our penalized regression models with different prior distributions predict the number of canyon heads in the hexagon tiles reasonably well with average prediction errors of RMSE = 1.2 and an explained variance of ∼50% (pseudo-R 2 of 0.5). Model performance, model parameters, and their posterior distributions remain largely unaffected by the choice of prior distribution ( Figures S5 and S6). Among the prior distributions, sampling with a lasso-shrinkage prior was most effective and thus we show the results of this model hereafter.
Predictions of our model are largely consistent with regional hotspots of shore-connected canyons along the Pacific coast of central South America and Central America, and the Mediterranean active margin (Figures 3a and 3b). Shore-connected canyons along the eastern Black Sea and some along the Indian-Ocean coast of the Arabic Peninsula are also identified. However, the model underestimates the frequency of shore-connected canyons along the Californian coast and the African passive margins (Figure 3). While the spatial patterns of modelled and actual shore-connected canyons are consistent, model residuals indicate high uncertainties in predicting spatial densities of canyons (Figure 3c). Figure 3d shows the posterior distributions of the regression coefficients of the 10 top-ranked predictors (Table S1). Variable ranking reveals the shelf gradient and the erodibility index as the top two predictors and the only predictors whose posterior distributions do not include zero in their 95% credible interval (regression coefficients of zero indicate no predictive value) (Figure 3d and Table S1). Two climatic parameters (Qw and TRMM 90/50 percentile) are ranked fifth and seventh and river bedload (QbBagnold) ranks tenth.

Data Limitations
Based on Bayesian reasoning and guarded against overfitting by applying penalized regression, we have identified credible predictors for shore-connected canyon occurrence along the world's coasts. However, we have not explored the fact that our predictors themselves and the response are prone to uncertainties. That our results are robust against the choice of prior distributions (and thus different degrees of uncertainty assigned to each parameter) suggests that our inferences are not strongly affected by these uncertainties.
Notwithstanding, there are additional uncertainties that cannot be quantified to date. For example, the onset of canyon incision has rarely been dated but can date back to several million years with repeated episodes of erosion and infilling (Maier et al., 2018;Mauffrey et al., 2017) and we have no global constraint on the age of canyon-head incision. Moreover, our predictors largely represent modern conditions, some of which may not represent active phases of canyon-head incision (e.g., Qs, rainfall, groundwater discharge, and wave height). Other variables, such as bedload, are difficult to determine, in particular on longer timescales (e.g., Nitsche et al., 2011). Additionally, canyons may incise along tectonic faults, reoccupy fluvial valleys on the shelf (Maier et al., 2018;Mauffrey et al., 2017), or preferentially incise along shelves built by erodible stratigraphy.
Moreover, longshore sediment supply can be crucial to canyon activity (e.g., Paull et al., 2013). While we included (water) depth of closure in the analysis, the width (and transport capacity) of the ocean-littoral cell is more appropriate to characterize longshore transport. This requires the combination of closure depth with BERNHARDT AND SCHWANGHART 10.1029/2020GL092234 10 of 15  high-resolution coastal bathymetry, data that is still unavailable globally (Bergsma et al., 2019). Sediment supply by glacial meltwater (Normandeau & Campbell, 2020) is also not included as we focus on latitudes 50°N-50°S.
Headward canyon erosion may be enhanced by wave-driven scouring which is determined by the orientation of wave crests with regard to the canyon (Smith et al., 2018). We included fresh SGD as a predictor, however, canyon formation may also be related to seepage of recirculated seawater (Pratson et al., 2007). Additionally, erosive sediment-gravity currents can be triggered by the seasonal downward flow of dense shelfwater (Canals et al., 2006), fostering shoreward canyon-head migration. Again, the above predictors were not included as they are unavailable for the spatial scale and extent of this study.
Finally, our analysis assumes that canyons are distributed independently from each other. However, canyon presence can influence fluid escape and thus the hydrology of the neighboring canyon (Orange et al., 1994). Such coupling mechanisms induce spatial autocorrelation (e.g., clustering) in canyon occurrence which ultimately may bias our modelling results. However, our Bayesian approach embraces the idea that uncertainties about model parameters depend on data availability and our posterior distributions can be updated once new data becomes available.

Controls on Shore-Connected Canyon Occurrence
Our findings support our first hypothesis that submarine canyons preferentially remain connected to the shoreline if shelves are steep and narrow (Figures 2b and 3d). Clearly, this situation minimizes the distance that a canyon head has to erode towards the shore over a single sea-level cycle. These findings, however, contrast those of Smith et al. (2017) along the US west coast where shelf width and gradient have no predictive power on canyon occurrence, a circumstance that is likely attributed to the low variability of shelf width (∼0 to <50 km) along this coast as compared to the global scale (∼0-493 km).
Shoreward canyon incision during sea-level rise occurs by mass wasting or erosive sediment-gravity flows. Steep shelves promote shoreward migration of canyon heads by accelerating flow velocities and bottom-shear stresses (e.g., Middleton, 1966), and hence, increased rates of canyon-floor downcutting. Downcutting, in turn, oversteepens canyon walls resulting in mass wasting and upslope-prograding failures (Densmore et al., 1997;Pratson and Coakley, 1996) or of canyon-thalweg knickpoints supporting shoreward canyon-head erosion (Guiastrennec-Faugas et al., 2020). Shelf-sediment failure in the canyon head can initiate turbidity currents, which further erode the canyon head (Pratson and Coakley, 1996). However, continental slope gradient has no credible influence on shore-connected canyon occurrence, suggesting a minor, if any, importance of retrogressive failure along steep continental slopes as a process in shore-connection maintenance.
Our global analysis partly supports the hypothesis of Smith et al. (2018) that shore-connected canyons occur preferentially offshore tectonically active regions underlain by resistant bedrock. Indeed, the three shore-connected canyon hotspots are located along tectonically active margins. However, topographic predictors indicating high onshore relief (catchment elevation and gradient) and uplift (channel steepness [ksn]) are ranked low in our model (Table S1) or even suggest a negative relationship (Figure 3d). Such relation is counterintuitive but may be explained by the effect of vast low-relief alluvial plains (especially in large catchments) that often separate tectonically active regions from the shore on the geomorphometric variables. The erodibility index holds rank 2 and only two shore-connected canyons occur at erodibility indices >2 (Figures 2b and S2a), which supports the hypothesis that resistant lithologies are important drivers on the incision of canyon heads. In fact, that shore-connected canyons are linked to catchments with high water discharge (rank 5) and availability of river bedload (ranks 10 and 11), suggests that river bedload may promote incision into the underlying shelf stratigraphy (Cook et al., 2013;Smith et al., 2017).
Our model underestimates the occurrence of shore-connected canyons along the Californian coast and misses some of the most prominent sediment conduits along passive margins (e.g., Congo and Capbreton canyons). While the detailed reasons for underestimation remain elusive, the misfit between model and observation along the Californian margin highlights the potential role of longshore sediment transport within littoral cells, an important process documented for several Californian canyons (Covault et al., 2011;Paull et al., 2005Paull et al., , 2013Romans et al., 2009), but which we cannot include in our predictors. The singular occurrence of canyons such as the Congo also remains unclear and we can only speculate that these features occur due to conditions that could not be considered in our global assessment, such as the reoccupation of shelf-incised fluvial channels and/or underlying faults.

Conclusions
Canyon heads close to the 120m-depth contour, the shoreline during the LCM, are globally abundant (n = 798, along major continents [islands excluded] between 50°N and 50°S). Yet, presently, there are only 183 shore-connected canyons, most of which belong to three spatial hotspots: the Mediterranean active margin and the Pacific coasts of central South and Central America. Whereas non-shore-connected canyons can host active turbidity currents (Kudrass et al., 1998;Zhong & Peng, 2021), the efficiency of land-to-ocean material transfer and burial should be higher for shore-connected canyons. For example, Mediterranean shore-connected canyons efficiently transport litter and contaminants onto the basin floor (Dominguez-Carrió et al., 2020;Pierdomenico et al., 2019). Hence, we expect major terrestrially derived litter concentrations and organic carbon burial on the ocean floor where spatial hotspots of shore-connected canyons coincide with litter-producing urban areas and/or areas with high primary production and erosion rates (Hilton & West, 2020).
Using Bayesian penalized regression, we predict the spatial patterns of these hotspots using a subset of predictors. We show that shore-connected canyons prevail along margins where shelves are narrow and steep and the onshore catchments expose resistant bedrock. Hence, on a global scale, low shelf width and high shelf gradients precondition the maintenance of canyon-head connectivity to the shore as the distance a canyon head has to erode shoreward during sea-level rise is minimal and high shelf gradients foster erosive sediment-gravity flows. Our analysis confirms previous studies (Smith et al., 2017(Smith et al., , 2018 which underscore the role of resistant bedrock exposed in onshore catchments in promoting submarine canyon incision. Together with high water discharge, these catchments deliver coarse-grained bedload which acts as a tool to erode the canyon head and floor.