Confinement width and inflow-to-sediment discharge ratio control the morphology and braiding intensity of submarine channels: insights from physical experiments and reduced-complexity models

Submarine channels conveying sediment gravity flows are often topographically confined, but the effect of confinement width on channel morphodynamics is incompletely understood. We use physical experiments and a reduced-complexity model to investigate the effects of confinement width (B) and the inflow-to-sediment discharge ratio ($Q_{\mathrm{in}}/Q_{\mathrm{s}}$) on the evolution of submarine braided channels. The results show that a larger confinement width results in increased active braiding intensity (BI_A) and that BI_A takes longer to stabilize (i.e., a longer critical time; t_c). At a fixed confinement width, a higher $Q_{\mathrm{in}}/Q_{\mathrm{s}}$ slightly decreases the BI_A. Digital elevation models of difference (DoD) of the experiments allow measurement of the morphological active width (W_a) of the submarine channels and the bulk morphological change (V_bulk) within an experiment, defined as the sum of total erosion and deposition. We find that W_a and V_bulk are proportional to B. We further confirm that BI_A is proportional to both dimensionless sediment–stream power (${\mathit{\omega }}^{**}$) and dimensionless stream power (ω^*). These trends are consistent for submarine braided channels both with and without confinement width effects. Furthermore, we built a reduced-complexity model (RCM) that can simulate flow bifurcation and confluence of submarine braided channels. The simulated flow distribution provides reliable predictions of flow depth and sediment transport rate in the experiments. Using kernel density estimation (KDE) to forecast the probability and trends of cross-sectional flow distribution and corresponding BI_A under extreme events, we find that skewness of the flow distribution decreases as discharge increases. The development of braided submarine channels, shown here to extend to conditions of topographic confinement, suggests that factors not modeled here (e.g., fine sediment) may be necessary to explain the abundance of single-thread submarine channels in nature.

1 Introduction

High-resolution bathymetry and seismic data have revealed unprecedentedly detailed channels on submarine fans and deep-sea plains (Deptuck et al., 2007; Babonneau et al., 2010; Janocko et al., 2013). These submarine channels are mainly formed by sediment gravity flows, and the corresponding channel morphology and stratigraphy resemble their fluvial counterparts (Lajeunesse et al., 2010; Peakall and Sumner, 2015; Jobe et al., 2016). Most past studies of submarine channels have focused on meandering channels (Imran et al., 1999; Peakall et al., 2000, 2007; Keevil et al., 2006; Deptuck et al., 2007; Straub et al., 2008; Sylvester et al., 2011). However, laboratory studies and field observations suggest that braided channels also develop in the submarine environment and may develop morphologies similar to braided rivers (Belderson et al., 1984; Hesse et al., 2001; Foreman et al., 2015; Lai et al., 2017; Limaye et al., 2018).

For instance, experiments of laterally unconfined submarine braided channels were first confirmed at two laboratory facilities (Foreman et al., 2015) and demonstrated that density currents readily produce submarine braided channels for flow aspect ratios (depth to width) similar to those of fluvial braided rivers. The authors confirmed that the stability theory of river planform morphology (Parker, 1976) successfully describes submarine braided channels in both experiments and the field. Subsequent experiments (Lai et al., 2017) further demonstrated that laterally unconfined submarine braided channels (Fig. 1a) have similar responses to fluvial braided rivers (Egozi and Ashmore, 2009). In this scenario, both the active braiding intensity (BI_A) and total braiding intensity (BI_T) increase with increasing inflow discharge and bed slope. BI_T is defined as the total number of channels in the cross section, while BI_A is a subset of the cross-sectional channel numbers, which reflects channels that can transport bed load, influence the channel morphology, and relate to channel bifurcations and avulsions (Ashmore, 2009; Egozi and Ashmore, 2009). Past fluvial and submarine braiding experiments (Ashmore, 1991; Egozi and Ashmore, 2009; Bertoldi et al., 2009; Lai et al., 2017) show that BI_A scales linearly with dimensionless stream power (ω^*). While the recently proposed entropic braiding index (Tejedor et al., 2022) accounts for the variation in the width of individual channel threads, in this study we use the traditional braiding index for ease of comparison with respect to these previous studies.

Figure 1 Submarine braided channels (a) in a laboratory-scale sedimentary basin, without lateral confinement (Lai et al., 2017), and (b) on the Baia di Levante submarine fan (Vulcano Island, Italy), with erosive banks as lateral confinement (modified from Romagnoli et al., 2012).

In natural rivers, lateral confinement reduces braiding intensity (Garcia Lugo et al., 2015; Carbonari et al., 2020). In submarine settings, canyon walls, erodible terraces, and basement ridges could also limit the development of submarine braided channels. For instance, the Orinoco deep-sea fan is the first reported field evidence of submarine braided channels and has a width-to-depth ratio of around 60 to 70 (Belderson et al., 1984; Ercilla et al., 1998). Further downstream of these submarine braided channels, confinement by the Barbados deformation front and a basement ridge forces the channels to transform into a single channel (Belderson et al., 1984; Callec et al., 2010). Along the Sicilian margin, submarine braided channels occur on the Stromboli slope valley (southeastern Tyrrhenian Sea) that is confined by a terrace and valley wall (Gamberi et al., 2011). In the same region, the Baia di Levante fan (Vulcano Island, Italy), a modern submarine volcaniclastic fan, was formed by the repetitive stacking of gravity flows (Romagnoli et al., 2012), with clear erosive banks confining the observed submarine braided channels (Fig. 1b).

No studies have investigated the evolution of submarine braided channels restricted by lateral confinement yet. In this study, we investigate how submarine braided channels respond to confinement widths (B) and the inflow-to-sediment discharge ratio ($Q_{\mathrm{in}}/Q_{\mathrm{s}}$) by using physical experiments and a reduced-complexity model. Our aims are (1) to obtain high-resolution digital elevation models (DEMs) for the evolving topography of submarine braided channels, (2) to analyze the effect of the confinement widths on the active braiding intensity, (3) to develop a reduced-complexity model for predicting discharge distribution patterns, and (4) to establish scaling relationships across fluvial and submarine braided channel systems.

2 Methods

2.1 Experiments

The experimental basin (1.8 m long, 0.55 m wide, and 0.05 m deep) was submerged into a water tank (2.3 m long, 1 m wide, and 0.65 m deep) to investigate the effect of different confinement widths on the evolution of submarine braided channels (Fig. 2). The water tank provided still, ambient water to simulate a deep-sea environment. Three drains were set at the downstream end of the tank to remove turbid water. The initial bed slope (S₀) of the basin is adjustable. Prior to the experiments, plastic sand (specific gravity = 1.5, median grain size d₅₀= 0.34 mm, and uniformity coefficient $C_{\mathrm{u}}=d_{\mathrm{60}}/d_{\mathrm{10}}=$ 1.64) was formed into a planar bed within the basin. Then a straight channel was excavated in the middle, with a fixed depth of 0.015 m and a prescribed confinement width. At the upstream end, the same dry plastic sand was provided steadily by a motor-controlled conveyor belt to simulate bed load sediment input. Dyed and saturated saline inflows (ρ_in= 1200 kg m⁻³) were supplied upstream to simulate long-lived, unconfined hyperpycnal flows or mud-rich turbidity currents (Métivier et al., 2005; Spinewine et al., 2009; Sequeiros et al., 2010; Weill et al., 2014; Foreman et al., 2015; Lai et al., 2016, 2017). Check dam partitions at the inflow exits allowed the saline discharges and dry plastic sands to be fully mixed before entering the basin.

Figure 2 (a) Experimental setup for studying the evolution of submarine braided channels with different confinement widths and inflow-to-sediment discharge ratios. (b) Inflow width, confinement width, and valley width of Runs A1 and B1. (c) Inflow width, confinement width, and valley width of Runs A2 and B2. (d) Inflow width, confinement width, and valley width of Runs A3 and B3.

Six runs of experiments were conducted to investigate the effects of a different confinement width (B) and inflow-to-sediment discharge ratio ($Q_{\mathrm{in}}/Q_{\mathrm{s}}$) on the evolution of the submarine braided channels. Combinations of confinement B (at 0.12, 0.24, and 0.48 m) and $Q_{\mathrm{in}}/Q_{\mathrm{s}}$ (at ∼ 60 and ∼ 90) encompass Series A and Series B (Table 1). The confinement widths were designed to be distinguished from previous experiments (Foreman et al., 2015; Lai et al., 2017) with a full valley width (W), and B in this study is based on the proportion of 21.8 %, 43.6 %, and 87.3 % of W. The inflow unit width discharge is calibrated by $q_{\mathrm{in}}=Q_{\mathrm{in}}/b$, where b is the inflow width (not confinement width B) because the submarine braided channels do not occupy the entire confinement width. The inflow-to-sediment discharge ratio ($Q_{\mathrm{in}}/Q_{\mathrm{s}}$) and inflow unit width discharge (q_in) were controlled at the same time. Based on the previously reported successful runs of Foreman et al. (2015) and Lai et al. (2017), $Q_{\mathrm{in}}/Q_{\mathrm{s}}=$ 60 and 90 are the two reasonable values of this ratio. The initial bed slope was kept at 0.114 (6.5^∘, which is the most effective bed slope for forming submarine braided channels reported by Lai et al., 2017). Each experimental run lasted 100 min and was divided into 10 successive stages (10 min for each stage). During the experiments, time-lapse images were obtained using a Nikon D7100 at intervals of 5 s. At the end of each stage, the inflow and sediment supply were stopped, and the topography was temporarily stabilized. Without draining the ambient water, the terrain surface was scanned. After that, a high-resolution (1 mm × 1 mm) digital elevation model (DEM), an orthorectified image, and a gradient map of each stage were constructed through a topographic imaging system (Lai et al., 2016, 2017).

Table 1 Summary of experimental conditions.

2.2 Automated measurement of channel position

Braiding intensity has been widely used in fluvial braided rivers (Egozi and Ashmore, 2008; Bertoldi et al., 2009; Egozi and Ashmore, 2009) and applied to submarine braided channels (Lai et al., 2017; Limaye et al., 2018). In past analyses of submarine braided channels (Lai et al., 2017), the braiding intensity required the manual interpretation of cross-sectional channel numbers. This approach is time-consuming and can be subjective. In this study, we developed an RGB (red, green, and blue) color separation algorithm to identify channel positions, calculate cross-sectional channel numbers, and obtain an objective, repeatable, and averaged braiding intensity for a selected color image (Fig. 3). For instance, a total of 1080 color orthophotos were analyzed using an algorithm based on Eq. (1) that obtains enhanced grayscale images.

$$\begin{array}{}\text{(1)}& \mathit{\alpha }\mathrm{G}-\mathit{\beta }\mathrm{R}-\mathit{\gamma }\mathrm{B},\end{array}$$

where α, β, and γ are the weighting coefficients to be calibrated by each run (e.g., α=1, β=0.8, and γ=0.35 for Run A2). R, G, and B represent the red, green, and blue bands of a color image, respectively. The weighting coefficients used to generate the grayscale images varied slightly between experimental runs due to different light settings and dye concentration. Then, the grayscale images were converted into binarized images by an adaptive threshold method (Bradley and Roth, 2007; Fig. 3b). After that, the cross-sectional illumination values (from x= 300 to 1500 mm, extracted every 1 mm) were smoothed by a Gaussian filter to locate the channel positions (Fig. 3c to Fig. 3e, showing only three example cross sections). After the above procedure, the channel positions and cross-sectional channel numbers could be automatically obtained (Fig. 3f). Then the average value would be taken as the representative active braiding intensity (BI_A) of that image. The results of the braiding intensity are presented in Sect. 3.2.

Figure 3 Procedure of image processing for identifying channel positions for each experiment. (a) Orthorectified color image. (b) Converted binarized image. (c–e) Extracted examples of the cross-sectional illumination values and captured channel positions. (f) Verification for the automatically captured channel positions and cross-sectional channel numbers.

2.3 Reduced-complexity model

Reduced-complexity numerical models are often used to simulate macroscale, long-term fluvial braided rivers (Murray and Paola, 1994, 1997; Thomas and Nicholas, 2002). The model proposed by Murray and Paola (1994, 1997) is the most influential work on reduced-complexity models. The model predicts the generic dynamics of stream braiding, including channel shifting, avulsion, and migration. Later research modified their flow routing schemes to make reduced-complexity models comparable to shallow water equation models (Thomas and Nicholas, 2002; Thomas et al., 2007; Nicholas, 2009; Williams et al., 2016). Recently, the CAESAR-LISFLOOD model (Coulthard et al., 2013) was applied to map experimental submarine braided channels and bar geometry (Limaye et al., 2018). This analysis showed that compared to a subaerial braided channel with the same water and sediment fluxes, submarine braided channels developed fewer and deeper channels near the inlet, with declining topographic variations occurring further downstream. In this study, we followed the numerical scheme of Thomas and Nicholas (2002) and developed a MATLAB-based reduced-complexity model, with modified hydraulic conditions for density currents. Our modifications assume large aspect ratios (width-to-depth ratio), reduced gravity and underflows that entrain negligible ambient water, and that the flows of submarine braided channels may resemble fluvial braided rivers (Foreman et al., 2015; Lai et al., 2017; Limaye et al., 2018).

The model starts calculations when the discharge is provided in the cell at the upstream boundary. Then the discharge is transferred downstream according to the calculation path of each row, and the model step ends when the calculation reaches the downstream boundary. The discharge is distributed from each upstream cell to the five immediate downstream cells (Thomas and Nicholas, 2002), rather than three cells used by Murray and Paola (1994, 1997). The discharge could be transferred laterally up to 60^∘ in the downstream direction. The allocated discharge in the five downstream cells is given by the following:

$$\begin{array}{}\text{(2)}& q_i={\displaystyle \frac{P_i}{{\sum }_{i=\mathrm{1}}^{\mathrm{5}}{P}_i}}{q}_{\mathrm{0}},\end{array}$$

where q₀ is the outflow discharge of the upstream cell. P_i is the routing potential, which determines the proportion of discharge that can be allocated to the five downstream cells. In the five downstream cells, the P_i of each cell was calculated from the upstream cell, starting from its water depth h₀. Assuming that the density current is uniform, the flow depth at steady state can be approximated by the following (Ippen and Harleman, 1952; Lofquist, 1960):

$$\begin{array}{}\text{(3)}& h_{\mathrm{0}}=C{\left({\displaystyle \frac{{\mathit{\mu }}_{\mathrm{in}}{q}_{\mathrm{0}}}{Jg\left({\mathit{\rho }}_{\mathrm{in}}-{\mathit{\rho }}_{\mathrm{a}}\right)S}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}},\end{array}$$

where C is a coefficient to be calibrated by each experiment (C= 1.2 for this study). ${\mathit{\mu }}_{\mathrm{in}}=\mathrm{1.68}\times {\mathrm{10}}^{-\mathrm{3}}$ kg m⁻¹ s⁻¹ is the dynamic viscosity of saline water at 20^∘, S is the bed slope, ρ_in is the inflow density, ρ_a is the density of ambient water, and J is a dimensionless parameter, which is proportional to the ratio of the interface velocity (V_i) to the maximum velocity (V_max) of the density current, which is given by the following:

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{V_i}{V_{\max }}}={\displaystyle \frac{\mathrm{12}J-\mathrm{1}}{\mathrm{12}{J}^{\mathrm{2}}+\mathrm{4}J+\frac{\mathrm{1}}{\mathrm{3}}}}.\end{array}$$

For a uniform density current at the steady state, $V_i/V_{\max }\cong$ 0.59 (suggested by Ippen and Harleman, 1952). Therefore, the calculated J value is 0.138 for this study. Flow depths (h_i) at each of the five downstream cells were approximated from the following:

$$\begin{array}{}\text{(5)}& h_i=h_{\mathrm{0}}+z_{\mathrm{0}}-z_i-S_{\mathrm{0}}dx,\end{array}$$

where z₀ is the bed elevation of the upstream cell, z_i is the bed elevation of each of the five downstream cells, and dx is the distance from the upstream cell to each of the five downstream cells. If the h_i calculated by Eq. (5) is a positive value, then the flow can be routed; if h_i is a negative value, then the flow cannot be routed. Therefore, the values of h_i affect the calculation method of P_i, which could be divided into the following two cases:

1. When at least one of the five downstream cells has h_i>0, then P_i can be calculated. In the remaining cells, if its h_i≤0, then its corresponding routing potential will be set to 0. P_i is calculated as follows:

   $$\begin{array}{}\text{(6)}& P_i=h_i^{\mathrm{1.67}}{S}_i^{\mathrm{0.25}},\end{array}$$

   where S_i is the local bed slope between the cells that flow is being routed from and to. The slope exponent of 0.25 was suggested by Thomas and Nicholas (2002).

2. If all five downstream cells have h_i ≤ 0, then the flow is assumed to be critical. This behavior often happens when density currents flood over sandbars, which is more likely to occur in submarine braided channels than in fluvial braided rivers (Sequeiros, 2012). Flow depths (h_i) at each of the five downstream cells are then approximated by the following:

   $$\begin{array}{}\text{(7)}& h_i=h_{\max }+z_{\min },\end{array}$$

   where h_max is the maximum depth, $h_{\max }={\left(q_{\mathrm{0}}/\sqrt{g^{\prime }}\right)}^{\mathrm{0.67}}$; $g^{\prime }=\left({\mathit{\rho }}_{\mathrm{in}}-{\mathit{\rho }}_{\mathrm{a}}\right)g/{\mathit{\rho }}_{\mathrm{in}}$ is the reduced gravity, and z_min is the lowest bed elevation in the downstream cell. The corresponding P_i is as follows:

   $$\begin{array}{}\text{(8)}& P_i=h_i^{\mathrm{1.5}}.\end{array}$$

Together, the above rules constitute our reduced-complexity model for simulating the flows of submarine braided channels. In our simulations, the initial bed condition is imported from the experimental DEM. The simulated inflow discharge is uniformly distributed at the upstream cells based on the experimental total inflow discharge. Note that our reduced-complexity model only simulates the underflows at the stable phase (Stage 6 to Stage 10); i.e., the morphological changes are excluded. Additionally, extra simulations are performed with double-inflow discharge, based on the experimental DEM at the final stage. These simulations are used for testing whether the linear relationship still holds between dimensionless stream power and active braiding intensity under extreme events. These results will be addressed in Sects. 3.4 and 4.

3 Results

3.1 Morphological evolution of submarine braided channels

Physical experiments document the initiation and evolution of submarine braided channels in response to set, known input parameters and conditions (Fig. 4). For instance, in Run A2, the initial confinement width was set at B= 0.24 m, and the inflow-to-sediment discharge ratio was $Q_{\mathrm{in}}/Q_{\mathrm{s}}=$ 58.7. By t= 1200 s, two small channels and a few sandbars appeared within the straight channel (Fig. 4b). At t= 2400 s, multiple channels and sandbars gradually developed downstream. At t= 3600 s, the flows showed frequent bifurcation and convergence. At t= 4800 s, slight bank erosion appeared on both sides, and some channels were abandoned. At t= 6000 s, channels frequently changed their routes, and sandbars grew and elongated. The DEM hillshades clearly document the channel relief and bar thickness of submarine braided channels (Fig. 5). These high-quality DEMs provided the basis for subsequent morphometric analyses.

Figure 4 Evolution of submarine braided channels of Run A2 from t= 0 to 6000 s. The confinement width is B= 240 mm. The inflow unit width discharge is q_in= 55.68 mm² s⁻¹. The inflow-to-sediment discharge ratio is $Q_{\mathrm{in}}/Q_{\mathrm{s}}=$ 59.

Figure 5 DEM hillshades for Run A2 from t= 0 to 6000 s.

Amongst Series A and Series B (Fig. 6), under the same $Q_{\mathrm{in}}/Q_{\mathrm{s}}$ condition, a higher confinement width produces more channels and a higher degree of braiding. On the contrary, when the confinement width is fixed, a larger $Q_{\mathrm{in}}/Q_{\mathrm{s}}$ generates fewer channels but wider main channels. For instance, in Series A ($Q_{\mathrm{in}}/Q_{\mathrm{s}}=$ 60), Run A3 (Fig. 6c) has more channels and a higher degree of braiding. Similarly, in Series B ($Q_{\mathrm{in}}/Q_{\mathrm{s}}=$ 90), Run B3 (Fig. 6f) is more braided. However, when the confinement width is fixed (e.g., B= 0.24 m; Fig. 6b and e), a larger $Q_{\mathrm{in}}/Q_{\mathrm{s}}$ would have fewer channels but wider main channels. These characteristics are also valid for confinements B= 0.12 m (Fig. 6a and d) and B= 0.48 m (Fig. 6c and f). Additionally, a larger $Q_{\mathrm{in}}/Q_{\mathrm{s}}$ also makes the sandbars more connected and facilitates small sandbars to merge into longer and larger sandbars (Figs. S1 to S10 in the Supplement; Movie S1 to Movie S6 provide more detailed evolution processes; Huang et al., 2023).

Figure 6 The last orthophoto at t= 6000 s for Series A (q_in≅ 55 mm² s⁻¹; $Q_{\mathrm{in}}/Q_{\mathrm{s}}\cong$ 60) and Series B (q_in≅ 84 mm² s⁻¹; $Q_{\mathrm{in}}/Q_{\mathrm{s}}\cong$ 90).

3.2 Evolution of active braiding intensity

The evolution of the active braiding intensity reflects the critical time for submarine braided channels to reach a stable phase and the influence of the confinement width and inflow-to-sediment discharge ratio (Fig. 7). For instance, the evolution of the braiding index over time for each run included two phases, namely (1) a rising phase and (2) a stable phase. We interpret the time for development of the stable phase as being a critical time (t_c). Our results showed that t_c is proportional to B. For example, when B= 0.12 m, t_c= 2400 s; when B= 0.24 m, t_c= 3600 s; and when B= 0.48 m, t_c= 4800 s. In addition, when $Q_{\mathrm{in}}/Q_{\mathrm{s}}$ is fixed, BI_A increases as B increases, e.g., Series A (Fig. 7a to c) and Series B (Fig. 7d to f). On the contrary, when B is fixed, BI_A at the stable phase decreases slightly as $Q_{\mathrm{in}}/Q_{\mathrm{s}}$ increases; e.g., BI_A values of Series B at t= 6000 s are all slightly smaller than that of Series A. Therefore, a larger confinement width postpones the critical time and increases the active braiding intensity at the stable phase. However, under the same confinement width, a larger inflow-to-sediment discharge ratio slightly decreases the BI_A.

Figure 7 Evolution of active braiding intensity (BI_A) for (a–c) Runs A1, A2, and A3 and (d–f) Runs B1, B2, and B3. t_c is the critical time for BI_A to reach a stable phase.

3.3 Volume change and morphological active width

The DEM of differences (DoDs) demonstrate the sediment distribution patterns of erosion and deposition for submarine braided channels (Fig. 8). The results showed that the confinement width and bank erosion were not proportional for Series A or Series B. For instance, when B= 0.24 m, the bank erosion on both sides was the most significant, especially for Run B2 (Fig. 8e). On the contrary, when the confinement width is fixed, a larger $Q_{\mathrm{in}}/Q_{\mathrm{s}}$ would cause a more severe bank erosion. Moreover, in Series B (Fig. 8d to f), the areas of erosion and deposition became more continuous and contiguous. In Series A (q_in≅ 55 mm² s⁻¹, $Q_{\mathrm{in}}/Q_{\mathrm{s}}\cong$ 60), most of the submarine braided channels develop within the given confinement width without much widening to accommodate the inflow. However, in Series B (q_in≅ 84 mm² s⁻¹, $Q_{\mathrm{in}}/Q_{\mathrm{s}}\cong$ 90), the submarine braided channels would widen the given confinement width to accommodate the larger inflow and cause stronger bank erosion. Besides, the DoDs also reveal the morphological active widths (W_a) and bulk changes (V_bulk) for each run. The relationship between W_a and V_bulk for both fluvial and submarine braided channels will be discussed in Sect. 4.

Figure 8 The DEM of difference (DoDs) for Series A (q_in≅ 55 mm² s⁻¹; $Q_{\mathrm{in}}/Q_{\mathrm{s}}\cong$ 60) and Series B (q_in≅ 84 mm² s⁻¹; $Q_{\mathrm{in}}/Q_{\mathrm{s}}\cong$ 90). The time duration for DoD calculation is dt= 4800–6000 s (last 20 min).

The normalized cross-sectional volume changes can further quantify the heterogeneity of the DoDs (Fig. 9). We extracted the DoD values on each cross section (from x= 0 to 1800 mm; every 1 mm), calculated the sum of its volume (positive value represents deposition; negative value indicates erosion), and normalized it by the absolute value of total volume change (the sum of the entire DoD values). The results show that the medians calculated by all runs are greater than 0, indicating that the DoD of each run is net depositional. When B= 0.12 m (Run A1 and Run B1), the DoDs are relatively homogeneous; i.e., the maximum and minimum values are closer to the median, showing a symmetrical distribution, but there are more outliers. When B= 0.48 m (Run A3 and Run B3), the DoDs are relatively heterogeneous; i.e., the maximum and minimum values are far from the median, showing an asymmetric distribution but fewer outliers. These trends hold for both $Q_{\mathrm{in}}/Q_{\mathrm{s}}=$ 60 (Series 1) and 90 (Series 2). Therefore, we propose that confinement width and heterogeneity of DoD are positively correlated.

Figure 9 Normalized cross-sectional volume changes for all runs.

3.4 Simulated flow discharge and sediment transport rate

The reduced-complexity model for simulating submarine braided channels captured the pattern of flow bifurcation and confluence (Fig. 10). Modeled flow occupied the main channels that developed in each run of the experiment. In Series A (Fig. 10a to c), a larger confinement width increased the number, width, and discharge of the main channels. In Series B (Fig. 10d to f), the main channels became wider and more significant. The simulated flow characteristics are similar to those qualitatively observed in the experiments. However, validating the reduced-complexity model by directly comparing the measured flow depths is unfeasible because the averaged depth of the experimental density currents is less than 3 mm. Therefore, we converted the simulated and experimental flow maps (Figs. S11 to S16 in the Supplement) into binarized images and compared their pixel differences. We examined the cross-sectional differences between the experimental and simulated flows every 1 mm along the x direction. The results (Fig. S17) showed that the errors are less than 30 % for B= 0.48 m (Run A3 and Run B3) and even less than 6.5 % for B= 0.12 m (Run A1 and Run B1). This indicates that our reduced-complexity model is a reliable and applicable tool for simulating submarine braided channels.

Figure 10 Submarine braided channels simulated by the reduced-complexity model for each run at t= 6000 s.

The simulated flows extracted at the middle reach illustrate the discharge, flow depth, and sediment flux distribution patterns affected by different confinement widths (Fig. 11). For instance, for Series A, the modality of the main channels belongs to a unimodal distribution (Fig. 11a to c); for Series B, the modality of main channels includes bimodal and multimodal distributions (Fig. 11d to f). Combining the simulated flow discharges and the bed load equation of Schoklitsch (1950), the sediment flux for each of the five downstream cells could be roughly estimated by the following:

$$\begin{array}{}\text{(9)}& q_{\mathrm{bi}}=\mathrm{2500}{S}_i^{\mathrm{3}/\mathrm{2}}\left(q_i-q_{\mathrm{c}}\right),\end{array}$$

where S_i is local bed slope, and q_i is the flow discharge (m² s⁻¹). q_c is the critical discharge (m² s⁻¹) and is given by the following:

$$\begin{array}{}\text{(10)}& q_{\mathrm{c}}=\mathrm{0.26}{\left(G-\mathrm{1}\right)}^{\mathrm{5}/\mathrm{3}}\left({\displaystyle \frac{d_{\mathrm{s}}^{\mathrm{3}/\mathrm{2}}}{S_i^{\mathrm{7}/\mathrm{6}}}},\right),\end{array}$$

where G is the specific gravity of the sediment (G= 1.5 for plastic sand), and d_s is the medium grain size (d_s= 0.34 mm for this study), which is suitable for a grain size between 0.31 to 7.02 mm. Note that there is no sediment flux when q_i < q_c.

Figure 11 Simulated discharge, flow depth, and sediment flux extracted at the middle reach (x= 900–1000 mm) for each run. The initial bed surface is based on the experimental DEM at the final stage.

Kernel density estimation (KDE) confirms that the probability distribution of simulated flows are robust and consistent. Under extreme events, the KDE provides the predicted pattern of the cross-sectional discharge and the corresponding BI_A (Fig. 12). For instance, we extracted the simulated flow discharge at the middle reach (x= 950 to 1050 mm) to perform probability density function (PDF) and KDE statistical analyses. We excluded the data for which its unit discharge is less than 50 mm² s⁻¹ to avoid too many outliers. The KDE shows that the PDFs of simulated discharges are all positively skewed and less sensitive to the confinement width. On the contrary, when doubling the inflow discharge, the skewness of KDE decreases, thus allowing the given discharges be allocated over a broader unit discharge. Therefore, the number of channels and corresponding BI_A all increased. Our reduced-complexity model appears consistent with the behavior of submarine braided channels, although the model still requires calibration and validation with more field bathymetric surveys.

Figure 12 Kernel density estimation (KDE) for the probability density function (PDF) of the simulated flow distribution for each run. Thick colored lines represent the KDEs of different runs. Dashed black lines are the KDEs of the double-inflow discharge.

4 Discussion

The purpose of this study is to understand the effects of the confinement width and inflow-to-sediment discharge ratio on the evolution of submarine braided channels using physical experiments and reduced-complexity models. We propose that confinement width has a stronger influence than the inflow-to-sediment discharge ratio on the morphological evolution of submarine braided channels. The scaling relationships established by experiments and simulations can help us better interpret the evolution of field-scale submarine braided channels.

First, we find that both active width (W_a) and bulk change (V_bulk) are proportional to the confinement width (B). Under a fixed reach length, a wider B would generate a larger active area. Then, the larger active area would give a larger active width. By definition (Ashmore et al., 2011; Peirce et al., 2018), the bulk change is the sum of the absolute values of the erosion volume (V_e) and deposition volume (V_d) in a DoD. A wider B increases the non-zero areas in a DoD, therefore increasing the V_bulk. For example, the W_a= 0.12–0.14 m for B= 120 mm (Fig. 13), the W_a= 0.24–0.28 m for B= 240 mm, and the W_a= 0.48–0.49 m for B= 480 mm. Our results show that W_a is proportional to V_bulk, which agrees with the response of the fluvial braided rivers (gray dots in Fig. 13). Even if the Reynolds number (Re) and densimetric Froude number (Fr_d) of our submarine braided channels indicate that flows in these channels are laminar and subcritical (Table 2), as long as the depth-to-width ratio ($h/B$) decreases and slope-to-densimetric Froude number ratio ($S/{\mathit{Fr}}_{\mathrm{d}}$) increases, submarine braided channels will naturally form (Fig. 14). Our results agree with the morphological classification of fluvial rivers (Parker, 1976). In short, the active width and bulk change in both submarine and fluvial braided channels share the same trend. This implies that the two braiding systems may have similar physical mechanisms.

Figure 13 The trend of active width (W_a) to bulk change (V_bulk). Colored dots represent the data of submarine braided channels (this study). Gray dots are the experimental data of the fluvial braided rivers (Peirce et al., 2018).

Table 2 Physical parameters and dimensionless numbers related to the submarine braided channels.

¹ Reynolds number, $\mathit{Re}=\frac{hu}{\mathit{\upsilon }}=\frac{Q}{\mathit{\nu }}$, where h is the flow depth of the saline underflow, estimated from the experiments (h≈ 2 mm). $u=Q_{\mathrm{in}}/\left(hB\right)$ is the cross-sectional averaged flow velocity. Q_in is inflow total discharge. B is the confinement width. ν is the kinematic viscosity of water (ν= 10⁻⁶ m² s⁻¹ at 20 ^∘C). ² Densimetric Froude number, ${\mathit{Fr}}_{\mathrm{d}}=\frac{u}{\sqrt{g^{\prime }h}}$, and $g^{\prime }=g({\mathit{\rho }}_{\mathrm{in}}-{\mathit{\rho }}_{\mathrm{a}})/{\mathit{\rho }}_{\mathrm{in}}$ is the reduced gravity. ρ_in is the density of inflow (ρ_in= 1200 kg m⁻³). ρ_a is the density of ambient water (ρ_a= 1000 kg m⁻³). ³ Dimensionless critical time, $t_{\mathrm{c}}^{*}=t_{\mathrm{c}}{Q}_{\mathrm{in}}/d_{\mathrm{s}}^{\mathrm{3}}$, where t_c is critical time. d_s is the sediment grain size (d_s= 0.34 mm). ⁴ Reduced bed shear stress, ${\mathit{\tau }}_{\mathrm{b}}=\left({\mathit{\rho }}_{\mathrm{in}}-{\mathit{\rho }}_{\mathrm{a}}\right)ghS$, where S is the bed slope. ⁵ Shields parameter, $\mathit{\theta }=\frac{{\mathit{\tau }}_{\mathrm{b}}}{({\mathit{\rho }}_{\mathrm{s}}-{\mathit{\rho }}_{\mathrm{in}})g{d}_{\mathrm{s}}}=\frac{\left({\mathit{\rho }}_{\mathrm{in}}-{\mathit{\rho }}_{\mathrm{a}}\right)ghS}{({\mathit{\rho }}_{\mathrm{s}}-{\mathit{\rho }}_{\mathrm{in}})g{d}_{\mathrm{s}}}=\frac{hS}{R^{\prime }{d}_{\mathrm{s}}}$, where ρ_s is the density of plastic sand (ρ_s= 1500 kg m⁻³). $R^{\prime }=({\mathit{\rho }}_{\mathrm{s}}-{\mathit{\rho }}_{\mathrm{in}})/\left({\mathit{\rho }}_{\mathrm{in}}-{\mathit{\rho }}_{\mathrm{a}}\right)$. ⁶ Boundary Reynolds number, ${\mathit{Re}}_p=\frac{u_{*}{d}_{\mathrm{s}}}{\mathit{\nu }}$, and $u_{*}=\sqrt{{\mathit{\tau }}_{\mathrm{b}}/({\mathit{\rho }}_{\mathrm{in}}-{\mathit{\rho }}_a)}$ is shear velocity. ⁷Dimensionless stream power, ${\mathit{\omega }}^{*}=\frac{\left({\mathit{\rho }}_{\mathrm{in}}-{\mathit{\rho }}_{\mathrm{a}}\right)QS}{{\mathit{\rho }}_{\mathrm{in}}{w}_{\mathrm{s}}{d}_{\mathrm{s}}^{\mathrm{2}}}$, and $w_{\mathrm{s}}=\frac{Rg{d}_{\mathrm{s}}^{\mathrm{2}}}{{\mathrm{C}}_{\mathrm{1}}\mathit{\nu }+(\mathrm{0.75}{\mathrm{C}}_{\mathrm{2}}Rg{d}_{\mathrm{s}}^{\mathrm{3}}{)}^{\mathrm{0.5}}}$ is the sediment settling velocity reported by Ferguson and Church (2004). $R=({\mathit{\rho }}_{\mathrm{s}}/{\mathit{\rho }}_{\mathrm{in}}-\mathrm{1}$) is the submerged relative density of the sediment (R= 0.25 in this study). C₁= 18 and C₂= 1 are two constants for typical natural sands. ⁸Dimensionless sediment–stream power, ${\mathit{\omega }}^{**}=\left(\frac{B}{h}\right)\left(\frac{Q_{\mathrm{s}}}{Q_{\mathrm{in}}}\right)S^{\mathrm{0.2}}$ (modified from Paola, 2001).

Figure 14 Morphological classification of fluvial rivers and submarine channels (modified from Parker, 1976, and Lai et al., 2017).

Second, we demonstrate that active braiding intensity (BI_A) is proportional to both the dimensionless sediment–stream power (${\mathit{\omega }}^{**}$) and the dimensionless stream power (ω^*), despite the dimensionless critical time ($t_{\mathrm{c}}^{*}$) of submarine and fluvial braided channels being different (Fig. 16). The dimensionless critical time is defined as $t_{\mathrm{c}}^{*}=t_{\mathrm{c}}{Q}_{\mathrm{in}}/d_{\mathrm{s}}^{\mathrm{3}}$ (Lai et al., 2017). We found that the $t_{\mathrm{c}}^{*}$ of lateral confined or unconfined submarine braided channels all fall around 10⁸ to 10⁹ (Table 2 and Fig. 16a), while the fluvial braided rivers fall around 10¹¹. We think that our smaller flume size, steeper bed slope, and lighter bed material together make the submarine braided channels reach $t_{\mathrm{c}}^{*}$ earlier than that of fluvial braided rivers. In addition, the Shields parameter (θ) corresponding to our sediment is about 0.46, which is much larger than the critical threshold for transport (∼0.04 to 0.06) and larger than what occurs in the fluvial braided rivers (Egozi and Ashmore, 2009; Fig. 15). The dimensionless sediment–stream power, defined as the combination of the width-to-depth ratio ($B/h$), sediment-to-inflow discharge ratio ($Q_{\mathrm{s}}/Q_{\mathrm{in}}$), and bed slope (S; Fig. 16b; modified from Paola, 2001), is proportional to the confinement width, sediment supply, and bed slope. When Q_s is fixed, a larger Q_in reduces ${\mathit{\omega }}^{**}$, making the corresponding BI_A become smaller. For instance, in Series B (q_in≅ 84 mm² s⁻¹; $Q_{\mathrm{in}}/Q_{\mathrm{s}}\cong$ 90), the submarine braided channels widen the imposed confinement width to accommodate the larger inflow, making the area of erosion and deposition in the DoDs more continuous and contiguous. This explains why increasing Q_in slightly reduces the BI_A of Series B at a stable phase. The ${\mathit{\omega }}^{**}$ is applicable to both unconfined and confined submarine braided channels (Fig. 16b). On the other hand, the dimensionless stream power (ω^*) considers the density differences in the underflow, grain size (d_s), and settling velocity (ω_s; Fig. 16c) but does not include the effects of B and Q_s. The ω^* works well for both unconfined fluvial braided rivers (Zanoni et al., 2007; Egozi and Ashmore, 2009; Bertoldi et al., 2009) and confined fluvial braided rivers (Garcia Lugo et al., 2015; Fig. 16c). In short, both ${\mathit{\omega }}^{**}$ and ω^* are proportional to BI_A. For an inverse problem, if BI_A is known, then these two relationships can be used to estimate the corresponding parameters (e.g., h, Q_s, Q_in, and density differences), which provides a way to apply this work to paleoflow reconstruction.

Figure 15 Shields diagram of fluvial rivers and submarine channels (modified from García, 2000).

Figure 16 (a) Dimensionless critical time ($t_{\mathrm{c}}^{*}$) of submarine and fluvial braided channels. (b) Relation of dimensionless sediment–stream power (${\mathit{\omega }}^{**}$) to active braiding intensity (BI_A) for two sets of submarine braiding experiments. (c) Relation of dimensionless stream power (ω^*) to active braiding intensity (BI_A) for two sets of submarine braiding experiments, four sets of fluvial braiding experiments, and one set of our reduced-complexity model.

Until now, the dimensionless parameters ($t_{\mathrm{c}}^{*}$, ${\mathit{\omega }}^{**}$, ω^*) calculated from experiments of submarine and fluvial braided channels have shown strong linear relationships and support our inference that the submarine braided channels on the Baia di Levante fan would continue to develop and prograde downstream. However, we think that the downstream lateral confinement may force the braided channels to transform into a single channel. The confinement width of these submarine braided channels gradually decreases from upstream to downstream and is well bounded by two erosive banks. The Baia di Levante fan was formed by the gradual stacking of gravity flows fed by repeated small-scale mass-wasting events rather than a single catastrophic collapse (Romagnoli et al., 2012). This modern field example agrees with the formation and morphological features observed in our laboratory experiments. Additionally, recent studies on the submarine flanks of Réunion Island also confirm that long-term erosive processes combined with a high sediment supply can form wide volcaniclastic deep-sea fans similar to the siliciclastic ones (Saint-Ange et al., 2011; Sisavath et al., 2011; Casalbore et al., 2021). Based on the reported bathymetry (Fig. 1b), when the water depth is 1000 m, the bed slope of the Baia di Levante fan is 4^∘, the confinement width is 3500 m (Romagnoli et al., 2012), and the number of submarine braided channels is about three. Assuming that these three submarine channels are active, we obtain BI_A= 3. Through the experimental scaling relationship (Fig. 16b), we obtain the corresponding dimensionless sediment–stream power as ${\mathit{\omega }}^{**}=$ 2.5. Therefore, we can estimate the reasonable flow depth range of the turbidity current by assuming $Q_{\mathrm{in}}/Q_{\mathrm{s}}$. For instance, when $Q_{\mathrm{in}}/Q_{\mathrm{s}}=$ 13, the estimated flow depth is 63.2 m. This depth agrees with the reported depth of the turbidity current (60 m) and the height of the erosive walls (70–80 m; Romagnoli et al., 2012). When $Q_{\mathrm{in}}/Q_{\mathrm{s}}=$ 9, the estimated flow depth is 91.4 m. It is possible to cause overspill and form sediment waves (or cyclic steps) outside the bank. Similar morphology appears on the left outer bank of the Baia di Levante submarine fan. Overall, the morphologic estimates derived from our proposed scaling relations are consistent with field observations, using plausible values for the inflow-to-sediment discharge ratio. Constraining these inputs for field cases would be major step toward predicting the morphology of submarine braided channels.

Additionally, the experimental results are intriguing in that they further emphasize the similarities in the propensity for and properties of braiding in fluvial and submarine systems (Lai et al., 2017). In isolation, the experimental results would predict that braided submarine channels are common in modern oceans; however, they are exceedingly rare (Foreman et al., 2015). Herein, we built upon previous experiments that documented the same relationship between dimensionless stream power and active braiding intensity among rivers and submarine channels (Lai et al., 2017) and expanded the experiments to show both the active width and bulk morphologic change increase as the confinement width (B) increases commensurate with experimental fluvial systems. Thus, braided submarine channels can be spontaneously generated both through self-channelization on an initially flat surface (Foreman et al., 2015; Lai et al., 2017; Limaye et al., 2018) and within a pre-existing confined channel (this study). Yet, in nature, both these bathymetric cases overwhelming result in sinuous, single-thread submarine channels. This implies there may be a missing attribute within the experiments that suppresses braiding in nature. Lai et al. (2017) suggested that the missing attribute was the enhanced deposition of suspended sediment loads within inactive channel threads during low-flow conditions. In nature, higher sediment concentrations within turbidity currents compared to river fluid flows likely facilitate greater aggradation and the in-filling of inactive channel threads. This process would limit the segmentation of the flow by bars and promote the development of a single thread over multiple threads. In the experiments, there is no suspended sediment load. This phenomenon may be exacerbated in field submarine channels for which the active width is confined by canyon or valley walls, compared to unconfined systems. In confined systems, the fine-grained turbidity currents that overspill channel margins cannot be advected across the overbank area to zones more distal from the channel. Instead, the fine-grained material is likely retained and forced to deposit within the active width. Interestingly, the Baia di Levante fan is comparatively coarse in its grain size distribution (Romagnoli et al., 2012) and, in the absence of significant silt and clay components, may represent a unique case in which a braided system may evolve. This suggests that provenance and the grain size distribution of sediment supply may play a nontrivial role in the submarine channel planform morphology. Finally, the observation that as the confinement width, B, increases, so the bulk morphologic change increase provides a testable hypothesis for stratigraphic studies of submarine channels. In other words, the increased bulk morphologic change should produce more heterogeneous channel sand body deposits that are particularly related to the abundance bar reactivation surfaces. It may be productive for future stratigraphic studies to quantify this abundance scaled to the channel geometry.

5 Conclusions

Physical experiments and a reduced-complexity model allow us to better explain the responses of submarine braided channels to confinement widths and inflow-to-sediment discharge ratios. The rules of stream power are applicable to field-scale submarine braided channels. Our main findings are as follows:

1. We confirm that a larger confinement width (B) postpones the critical time (t_c) and increases the active braiding intensity (BI_A) at the stable phase. At a fixed confinement width, a larger inflow-to-sediment discharge ratio ($Q_{\mathrm{in}}/Q_{\mathrm{s}}$) slightly decreases the BI_A.

2. The active width (W_a) and bulk change (V_bulk) are proportional to the confinement widths for submarine braided channels, which is in agreement with fluvial braided rivers.

3. The experimental and modeling results agree that active braiding intensity is proportional to both dimensionless sediment–stream power (${\mathit{\omega }}^{**}$) and dimensionless stream power (ω^*).

4. Kernel density estimation (KDE) helps us to visualize and predict the pattern of cross-sectional discharge distribution.

5. We infer that the submarine braided channels on the Baia di Levante fan may continue to develop and prograde downstream but that the downstream lateral confinement will eventually force the submarine braided channels to transform into a single channel.

6. Overall, the relations we develop from the experimental and modeling results are able to predict the first-order general morphological and sedimentological patterns for field-scale submarine fans.