Influence of alluvial slope on avulsion in river deltas

Changing hydrological regimes, sea-level rise, and accelerated subsidence are all putting river deltas at risk across the globe. One mechanism by which deltas may respond to these stressors is that of avulsion. Decades of delta avulsion studies have resulted in conflicting hypotheses as to whether avulsion timing and location are primarily controlled by upstream (water and sediment discharge) or downstream (backwater and sea-level rise) drivers. Here we use Delft3D morphodynamic simulations to test the upstream-influence hypothesis by varying the initial alluvial slopes upstream of a self-formed delta plain within a range ($\mathrm{1.13}\times {\mathrm{10}}^{-\mathrm{4}}$ to $\mathrm{3.04}\times {\mathrm{10}}^{-\mathrm{3}}$ m m⁻¹) that is representative of global deltas, while leaving all other parameters constant. Avulsion timing and location were recorded in six scenarios modelled over a 400-year period. We measured independent morphometric variables including avulsion length, delta lobe width, bankfull depth, channel width at avulsion, delta topset slope, and sediment load and compare these to natural and laboratory deltas. We find that larger deltas take more time to avulse, as avulsion timing scales with avulsion length, delta lobe width, and bankfull depth. More importantly, we find strong negative correlations between sediment load avulsion timescale and sediment load initial alluvial slope. Sediment load is directly dependent on the upstream alluvial slope, and increases in this slope raise transport capacity and introduce more sediment into a delta plain, leading to higher aggradation rates and, consequently, more frequent avulsions. These results introduce further debate over the role of downstream controls on delta avulsion.

1 Introduction

River deltas are home to ∼ 339 million people worldwide, are hotspots for biodiversity, and are crucial carbon sinks (Ericson et al., 2006; Hackney et al., 2020; Loucks, 2019; Shields et al., 2017; Syvitski and Saito, 2007). However, the geomorphic dynamism of river deltas has been, and continues to be, altered by amplifying stressors such as changing hydrological regimes, sea-level rise, and accelerated subsidence, putting human and other systems that rely on river deltas at considerable risk (Giosan et al., 2014; Syvitski et al., 2009; Tessler et al., 2015; Wallace et al., 2014). A frequently observed and geologically rapid mechanism by which deltas respond to these stressors is by river flow avulsing from one distributary channel into another. Delta avulsion location correlates with backwater length, slope break, and valley exit location measured from the shoreline (Ganti et al., 2016a; Hartley et al., 2017; Prasojo et al., 2022). Many studies have proposed different hypotheses for the main controls of delta avulsion frequency (e.g. Aslan et al., 2005; Brooke et al., 2020; Edmonds et al., 2009; Kleinhans and Hardy, 2013; Nijhuis et al., 2015; Slingerland and Smith, 2004), but there is currently no consensus over the conditions under which the various driving factors control this frequency.

During avulsion, flow is abruptly diverted out of an established river channel into a new course on the adjacent floodplain or delta plain (Jones and Schumm, 2009; Slingerland and Smith, 2004). When a delta channel avulses, the population, infrastructure, and economic activities on the delta plain may be at risk. Delta-lobe-scale avulsions may be considered rare, but this is partly due to anthropogenic controls on many delta channels preventing avulsion (e.g. built riverbanks), and unmodified systems can exhibit avulsion over decadal or shorter timescales, for example, once every 12 years in the Yellow River Delta (Jerolmack, 2009) or once every 4 years in the Sulengguole River, China (Li et al., 2022). Avulsions may be full, where the flow following a new course completely abandons its parent channel, or partial, in which only a portion of the flow is diverted (Slingerland and Smith, 2004). Avulsion is often effectively instantaneous but may also be gradual, as in the Rhine–Meuse delta, where one documented avulsion event took 1250 years to complete (Stouthamer and Berendsen, 2001). There are also several styles of avulsion: annexation, in which a pre-existing channel is reoccupied; incision, where a new channel is scoured into the floodplain surface as a direct result of the avulsion; and progradation, where extensive sediment deposition, such as a mouth bar, causes flow bifurcation and formation of a multi-channelled distributive network (Slingerland and Smith, 2004).

River deltas are initiated through repeated mouth bar deposition due to sudden expansion and deceleration of a sediment-laden jet of water entering relatively still water, usually a sea or lake (Bates, 1953; Edmonds et al., 2011; Kleinhans et al., 2013; Wright, 1977). Mouth bars grow in both upstream and downstream directions from the point of initiation. Once a mouth bar's aggradation reaches 40 %–60 % of the initial flow depth, it will stop growing because the sediment is advected around the bar rather than accelerated over it (Edmonds and Slingerland, 2007; Fagherazzi et al., 2015; Kleinhans et al., 2013). This cessation of bar growth is the point where avulsion by progradation or bifurcation starts in a river delta. Simultaneously, avulsion by incision takes place in the proximal parts of a delta plain when mouth bar deposition and stagnation induce parent channel backfilling or in-channel aggradation, triggering an avulsion to create a smaller distributive channel network by breaching the channel levee (Ganti et al., 2016a). The location of the point furthest upstream where a delta channel starts to avulse correlates with the location of a break (i.e. decrease) in bed slope (Prasojo et al., 2022; Ratliff et al., 2021), the limit of the backwater zone (Brooke et al., 2022; Ganti et al., 2016a), and the exit point from the river valley (Hartley et al., 2017).

A study of 105 global river deltas found that the strongest correlation is between the locations of avulsion nodes and the locations of the decreases in slope at the transition from an alluvial to a delta slope (Prasojo et al., 2022). Consequently, it is hypothesized that the slope of the alluvial river upstream of a delta controls the frequency of avulsion on delta plains, with steeper alluvial slopes leading to more frequent avulsions. This control is due to greater sediment transport capacity on steeper slopes (Bagnold, 1966) such that, subject to sediment availability, more sediment per unit width will be delivered to a delta plain where alluvial slopes are steeper. Assuming constant channel width and no subsidence, any reduction in stream power across the delta plain leads to aggradation, the rate of which will be greater when upstream sediment supply is higher, which in turn leads to increased avulsion frequency (Jerolmack and Mohrig, 2007; Mohrig et al., 2000). Alternatively, lower alluvial slopes are associated with lower sediment input flux and hence less frequent avulsion.

To test if the alluvial slope upstream of a delta controls the timing of avulsion on delta plains, we use Delft3D morphodynamic simulation software to (1) assess the effect of varying alluvial slopes upstream of a delta slope break on the avulsion timescale and to (2) investigate the primary controls over delta avulsion. Morphometric variables (delta lobe width, channel width at avulsion, avulsion length, topset slope, bankfull depth, and sediment supply) were measured at every time step during delta growth. These morphometric properties are measured as independent variables expected to covary with avulsion timescales. This investigation aims to (1) identify the role of alluvial slope upstream of delta plains on avulsion timescales, (2) explain the mechanisms by which the controlling variables determine avulsion timescale, and (3) compare avulsion timescales from this numerical model with an analytical solution and also with observations from natural and physical experimental river deltas. A robust understanding of these processes has practical implications due to their direct impact on coastal and inland flood risk in highly populated river deltas and contributes to a fundamental understanding of natural delta-building processes.

2 Methods

We designed a set of numerical experiments to model a natural-scale river delta (7.5×7.5 km, 300 by 300 computational cells, each 625 m²) using Delft3D (v.4.04.02) software. For comparability with previous studies, we adopted physical parameters used in similar Delft3D river delta models by Edmonds and Slingerland (2010) and Caldwell and Edmonds (2014). Model bathymetry was designed to accommodate the six alluvial slopes defined below as our model scenarios.

2.1 Scenario definition

The model used a range of alluvial slopes upstream of the delta's slope break (S_alluvial) (Fig. 1a), which are considered to be representative of natural river deltas (Fig. 1b). Representative percentiles of the ratio between S_alluvial and downstream delta topset slope (S_topset) were determined from 105 global river deltas measured by Prasojo et al. (2022) (Fig. 1b; Table 1). Percentiles of the $S_{\mathrm{alluvial}}/S_{\mathrm{topset}}$ ratio of 2.5, 10, 25, 50, 71, and 75 were used to define model scenarios. Model alluvial slopes were calculated from these ratios using a constant downstream slope (S_downstream=0.000375) similar to that of the Atchafalaya Basin, Mississippi Delta, Louisiana (Edmonds and Slingerland, 2010). During the simulation, both alluvial (S_alluvial) and downstream (S_downstream) slopes evolved and created delta plains with varying topset slope values (S_topset), as aggradation occurred (Video S1 in the Supplement). We defined equilibrium to have been reached when there was constant sediment discharge and channel depth at the model inlet after ∼ 3–6 d of simulation time.

Figure 1 (a) Schematic diagram of the model design. The alluvial slope of each run was calculated from six percentiles from the alluvial slope–topset slope ratios of modern river deltas shown in Fig. 1b. Initial downstream slope is kept constant at 0.000375, adopted from the downstream slope of the modern Mississippi Delta (Edmonds and Slingerland, 2010). (b) Distribution of the ratio between alluvial (S_alluvial) and topset (S_topset) slopes from 105 modern river deltas distributed across five climate regions. Ratios used for numerical model runs are indicated by dashed vertical lines. (c) Plan view of the model design. L_s and L_a are slope break and avulsion lengths, respectively. The non-erodible bed at 5 m above sea level represents non-erodible bedrock. (d) Schematic diagram of a river delta showing avulsion location, inlet sediment supply (Q_s), lobe width of each avulsion (B), avulsion length (L_a), and channel widths measured at avulsion (B_c), modified from Chadwick et al. (2020). Numbers near the shoreline represent the number of delta lobes that were used to measure B; e.g. B₄ in panel (d) represents the width of the fourth lobe built. (e) Schematic cross-section showing basin depth (H_b) and topset slope (S_topset). Parameters shown in Fig. 1d–e are measured at each time step during delta growth.

2.2 Model setup

We used Delft3D software to model six scenarios. Delft3D is a physics-based model that simulates hydro-morphodynamics and has been validated for a wide range of environments, including reach-scale meandering and braided rivers, estuaries, and self-formed river deltas (Edmonds and Slingerland, 2007, 2008; Geleynse et al., 2011; Leuven et al., 2023; Morgan et al., 2020; Nienhuis et al., 2018a; Nijhuis et al., 2015; Rossi et al., 2016; Williams et al., 2016). Delft3D solves the shallow-water equation and is integrated with the D-Morphology sediment transport and morphology module (Deltares, 2021). Delft3D calculates the flow velocity, sediment entrainment, transport, and deposition, and it updates bed levels at each computational time step by calculating bed sediment mass for each cell as the results of suspended and bedload transport divergence (Deltares, 2021).

Table 1 Numerical modelling scenarios as defined in Fig. 1.

The alluvial slope was the only controlled variable in our study, with all other setup and parameters the same as “scenario o” of Edmonds and Slingerland (2010) and Caldwell and Edmonds (2014) (Fig. 1c). The model domain was rectangular with four boundaries, the incoming river discharge being located at the “South” boundary of the model and the other three boundaries set to 0 m elevation above sea level (Fig. 1c). The constant incoming river discharge, set at 1050 m³ s⁻¹, was uniformly distributed across the 250 m wide inlet channel, and inlet sediment discharge was in equilibrium with transport capacity. Various alluvial slopes were achieved by changing the inlet channel elevation in each run while maintaining the receiving basin's depth. Consequently, sediment discharge varied in each run because of the varied alluvial slope as the main controlled variable in the experiments. Our modelled deltas closely represented natural deltas because the discharge ratio and the differences in bed heights between bifurcating distributary channels follow ranges similar to those reported for natural deltas (Edmonds and Slingerland, 2010). Sea level remained constant within the model, and no subsidence, tide, or wave effects were considered.

The model domain was 7.5 km × 7.5 km to avoid the delta plain extending across the model boundaries. Constant sediment grain size distributions were used throughout the model. Fine sand was introduced as non-cohesive sediment (D₅₀=125 µm with a Gaussian distribution), and medium-grained silt was introduced as cohesive sediment (D₅₀=30 µm). The critical bed shear stress for erosion was 0.10 N m⁻², and the model initially contained 5 m of fully mixed sediments. We introduced a slope break 1 km from the inlet boundary to drive delta formation in the model's initial bathymetry. Using the slope break and avulsion length scaling identified from measured global river deltas (Prasojo et al., 2022), the expected first avulsion node location should emerge in each scenario at around 2.2 km from the inlet (Fig. 1c). Other physical and numerical parameters were held constant across all scenarios (Table 2).

For 18 d simulation, the model produced one output every 480 min, resulting in 52 visualization outputs (i.e. maps) at the end of the simulations. Using a morphological scale factor (morfac) of 175, these 52 maps represented 3150 d (8.6 years) of prototype time. Because bankfull discharge occurs for ∼ 2 % of the time on average (Dunne and Leopold, 1978), 18 d of simulation thus represent around 430 years of “real” time (i.e. 8.6 years divided by 0.02).

Table 2 User-defined model parameters (adopted from Edmonds and Slingerland, 2010; Caldwell and Edmonds, 2014).

2.3 Surface morphological metrics

The model reached equilibrium after ∼ 3–6 d of simulation time when we began measurements of morphometric variables and avulsion timescale. Avulsion timescale, defined as the time in between successive avulsions, was empirically measured throughout all experiments after equilibrium was attained. Avulsions were defined when a distributary channel produced during delta formation changed its course and commenced deposition of a new delta lobe. We only considered avulsions caused by progradation or incision that are common in river deltas (Slingerland and Smith, 2004). The timing of each avulsion in the model was noted and converted to a “real” time as T_{a empirical}.

Numerous morphological surface metrics can be used to describe delta form. The surface metrics used here followed those used in an analytical solution for avulsion timescale (Eq. 4 from Chadwick et al., 2020), which utilized delta lobe width (B), channel width at avulsion (B_c), avulsion length (L_a), basin depth (H_b), magnitude of relative sea-level rise (RSLR; z), topset slope (S_topset), bankfull depth (h_c), and sediment supply (Q_s). Avulsion length, delta lobe width, channel width at avulsion, and delta topset slope were measured on all maps after equilibrium was reached. The delta lobe width, channel width at each avulsion node, and avulsion length were measured in QGIS from the georeferenced images produced by Delft3D (Fig. 1d, Table S1). Delta lobe width was measured as the maximum width of each lobe, while avulsion length was measured along the longest channel from the shoreline to the most upstream avulsion node every time new avulsion occurred in our models. Topset slope (S_topset) was calculated as the average slope of a delta plain. To measure S_topset, digital elevation models (DEMs) for each time step were first extracted from Delft3D and then cropped to only cover the delta plain. We then filtered the DEM to only include elevation, z between 0–5 m ($\mathrm{0}<z<\mathrm{5}$), to cover the part of the delta plain that is exposed above sea level (i.e. z=0 m) but below the non-erodible bed (i.e. z=5 m; Fig. 1d). After cropping and filtering, the DEM was then transformed to a slope raster defined as the change in elevation for each DEM cell in x and y directions with $S_{\mathrm{topset}}=\sqrt{{\left(\frac{\mathrm{d}z}{\mathrm{d}x}\right)}^{\mathrm{2}}+{\left(\frac{\mathrm{d}z}{\mathrm{d}y}\right)}^{\mathrm{2}}}$. The mean topset slope for each time step was then extracted from the slope raster as the slope values for each scenario showed unimodal distributions (Table S1; Fig. S1). Lastly, the sediment load (Q_s) at the channel inlet was obtained from the QUICKPLOT (v2.60.65942) Delft3D visualization software.

Bankfull depth (h_c) was calculated using Eq. (1) (Parker et al., 2007):

$$\begin{array}{}\text{(1)}& h_{\mathrm{c}}={\left({\displaystyle \frac{C_{\mathrm{f}}{Q}^{\mathrm{2}}}{g{B}_{\mathrm{c}}^{\mathrm{2}}{S}_{\mathrm{topset}}}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}},\end{array}$$

where C_f is a bed friction coefficient (–) equal to 0.002 for large lowland rivers (Parker et al., 2007), Q is bankfull discharge (m³ s⁻¹) equal to 1050 m³ s⁻¹, and g is gravitational acceleration (m s⁻²) equal to 9.81 m s⁻².

The avulsion timescale was calculated between each successive pair of avulsions observed in the model (T_{a empirical}) and was correlated with all the measured morphometric variables (e.g. Q_s, L_a, B_c, B, S_topset, S_alluvial, and h_c) from all post-equilibrium maps. Scatter plots and Pearson correlation coefficients (r) were used to assess the shape of relationships and potential dependencies between all variables.

2.4 Comparing numerical models to analytical solutions and natural deltas

We used the mass-balance-based analytical solution of Chadwick et al. (2020) to calculate expected avulsion timescales for our model conditions (Eqs. 2–5, Table S1). Measured independent surface morphological metrics were used in Eqs. (2)–(5) to calculate avulsion frequency (f_a) and timescale (T_a).

$$\begin{array}{}\text{(2)}& {\displaystyle }\begin{array}{rl}& f_{\mathrm{a}}={{\displaystyle \frac{\mathrm{1}}{T}}}_{\mathrm{a}}={\displaystyle \frac{\mathrm{1}}{\left(\mathrm{1}-{\mathit{\lambda }}_{\mathrm{p}}\right)}}\\ & {\displaystyle \frac{Q_{\mathrm{s}}}{\left(L_{\mathrm{a}}-D\right)\mathrm{BH}+\mathrm{DB}\left(H_{\mathrm{b}}+z+\frac{{\mathrm{DS}}_{\mathrm{topset}}}{\mathrm{2}}\right)}}\mathrm{if}D\ge \mathrm{0},\end{array}\text{(3)}& {\displaystyle }{f}_{\mathrm{a}}={\displaystyle \frac{\mathrm{1}}{T_{\mathrm{a}}}}={\displaystyle \frac{\mathrm{1}}{\left(\mathrm{1}-{\mathit{\lambda }}_{\mathrm{p}}\right)}}{\displaystyle \frac{Q_{\mathrm{s}}}{L_{\mathrm{a}}\mathrm{BH}}}\mathrm{if}D<\mathrm{0},\text{(4)}& {\displaystyle }D=(H-z)/S_{\mathrm{topset}},\text{(5)}& {\displaystyle }H=H^{\ast }{h}_{\mathrm{c}},\end{array}$$

where f_a is the avulsion frequency (year⁻¹), Q_s is the sediment load (m³ s⁻¹), λ_p is the sediment porosity (–), L_a is the avulsion length (m), D is the delta lobe progradation distance (km), B is the delta lobe width of each avulsion (m), H is the aggradation thickness necessary for avulsion (m), H_b is the basin depth (m), z is the magnitude of sea-level rise during an interavulsion period (m), S_topset is the topset slope (–), H^* is the avulsion threshold (–), and h_c is the bankfull depth (m) calculated using Eq. (1).

In calculating these analytical avulsion timescales, sensitivity analyses were undertaken using avulsion thresholds (H^*) of 0.2, 0.5, and 1.4, which are realistic for lowland deltas (Ganti et al., 2019), and D>0, since there is no allogenic forcing that would make the delta regress. The analytical avulsion timescales for $H^{*}=\mathrm{0.2}$, 0.5, and 1.4 are $T_{\mathrm{a}{H}^{*}=\mathrm{0.2}}$, $T_{\mathrm{a}{H}^{*}=\mathrm{0.5}}$, and $T_{\mathrm{a}{H}^{\ast }=\mathrm{1.4}}$, respectively (Table S1). Since sea level is constant in this study, the sea-level rise (z) is 0. Sediment porosity (λ_p) is assumed to be 0.4 (Jerolmack, 2009; Paola et al., 2011), the bed friction coefficient (C_f) is 0.002 for lowland rivers (Parker et al., 2007), and the constant bankfull discharge (Q) is 1050 m³ s⁻¹. Analytical avulsion timescales were then compared to avulsion timescales observed from 19 natural river deltas, 2 fan deltas, and 1 downscaled physical laboratory fan delta documented in Chadwick et al. (2020) and Jerolmack and Mohrig (2007), using topset slope values from Prasojo et al. (2022) (Table S2).

3 Results

Figure 2 shows the morphology of the deltas in each scenario at the final time step. Overall, the different alluvial slopes are associated with delta plains that exhibit different shoreline configurations, different numbers of active distributary channels, and slightly different delta plain sizes. One delta plain reached the model boundary (Run f, Fig. S2), and this scenario was repeated with a larger domain size, and the avulsion timescales and morphological metrics were observed from this larger domain.

Figure 2 (a–f) River deltas for each run at the final simulation time step. Run f uses a larger (12.5×12.5 km) model domain to avoid the delta plain reaching the model's boundary. Morphometric measurements for Run f were made on this larger domain.

As an example of avulsion in the model, Fig. 3 shows one avulsion that was initiated by an increase in in-channel deposition in a distributary channel around the avulsion node (Fig. 3a–g). In-channel aggradation sets up an avulsion to occur by increasing the likelihood of overbank flows because it elevates the water level and distributary channel bed above the surrounding floodplain (Fig. 3g). Aggradation continues until the levee is breached through erosion because the surrounding delta plain provides a lower path for a distributary channel to flow (Fig. 3b). This newly avulsed channel then distributes more water and sediment away from the initial distributary channel path (Fig. 3c–f). Note also that bifurcation at the delta shoreline occurred at the same time as the avulsion developed ∼ 1 km upstream (Fig. 3a). Bifurcation is initiated by a mouth bar deposition in a more distal part of the delta plain (Fig. 3a). The distributary channel feeding this mouth bar bifurcates once the depth of the mouth bar is ∼ 40 % of the initial basin depth, consistent with findings from Edmonds and Slingerland (2007). However, we are uncertain if the mouth bar deposition forces an upstream wave of in-channel deposition that leads to the channel being unstable enough to trigger an avulsion (Fig. 3a–g).

Avulsion occurs as often as every 8.27 years in all scenarios, with the maximum time between successive avulsions being 148.9 years (Run d). At the upper limit, avulsions occurred 14 times (runs a and f), and they occurred 6 times at the lower limits (Run c). A right-skewed distribution of avulsion interval was observed in most scenarios (Fig. S3), with medians from 8.3–16.5 years and means from 14.8–33.1 years. The median timing of avulsions remains unchanged over runs a–f (Kruskal–Wallis, p>0.05) (Fig. S3).

Figure 4 shows Pearson's correlations between observed avulsion timescales in the model (T_{a empirical}) and the independent morphometric variables measured in each time step. T_{a empirical} has a weak negative correlation with sediment load, $Q_{\mathrm{s}}(r=-\mathrm{0.06})$. The initial alluvial slope (S_alluvial; Table 1, Fig. 1), defined independently in our approach, controls the sediment load (Q_s; r=0.86) and is responsible in defining the delta avulsion length, L_a (r=0.58; Fig. 4), suggesting that higher transport capacity on a steeper alluvial slope may produce a longer (i.e. more protruding) delta, consistent with the findings from global river deltas (Prasojo et al., 2022). T_{a empirical} is also weakly and negatively correlated with bankfull depth, h_c ($r=-\mathrm{0.10}$), suggesting that bankfull depth may not be a good indicator of avulsion timescale. Sediment load (Q_s) influences the delta topset slope (S_topset; r=0.37) so higher riverine sediment loads produce steeper delta plains when the basin configuration is similar. Avulsion timescale (T_{a empirical}) shows no correlation with the delta size represented as delta avulsion length, L_a ($r=-\mathrm{0.01}$); delta lobe width, B ($r=-\mathrm{0.09}$); and channel width at avulsion node, B_c (r=0.04). Note that there are high degrees of autocorrelation within the morphometric parameters used here, for example, between bankfull depth (h_c) with channel width at avulsion, B_c ($r=-\mathrm{0.95}$), and with topset slope, S_topset ($r=-\mathrm{0.71}$).

Figure 3 Time series images from Run a showing the details of an avulsion process (a–f). Cross-section showing the evolution of bed and water level during the avulsion process (g).

Figure 4 Pearson correlations (upper-right panels) between avulsion timescale (T_{a empirical}) and independent morphometric variables, along with their distributions (diagonal panels) and correlations (lower-left panels). Units in this figure are years for T_{a empirical}; cubic metres per second (m³ s⁻¹) for Q_s; and metres (m) for B_c, h_c, and L_a. S_topset and S_alluvial are dimensionless. Note that h_c is autocorrelated with S_topset and B_c as shown in Eq. (1) and that S_alluvial is the initial alluvial slope available in Table 1 as the independent variable used to define the six experimental scenarios. Red lines on the correlation plots are LOWESS curves. Red lines on the histograms aid visualizing the distributions. Also note that T_{a empirical}, Q_s, S_topset, and L_a are on log scales.

As Fig. 4 groups the results together, Fig. 5 shows ordinary least-squares regressions between sediment load (Q_s), avulsion timescale (T_{a empirical}), initial alluvial slope (S_alluvial), and bankfull depth (h_c) produced from each scenario. Avulsion timescale is inversely and non-linearly correlated with sediment load (see R² values for each scenario and overall dataset at the bottom of Fig. 5), showing the influence of sediment load in defining avulsion timescale observed in our model (Fig. 5a). The higher the riverine sediment load brought into a delta plain, the faster avulsions occur, as shown in Fig. 5a. Conversely, the riverine sediment load is positively correlated with initial alluvial slope (Fig. 5b), showing how initial alluvial slope upstream of a delta plain defines how much riverine sediment is being transported to the delta plain. As we maintain the channel width upstream of a delta plain to be constant (Fig. 1c), a higher transport capacity in a steeper alluvial slope is able to bring more sediment load (Fig. 5b). Lastly, Fig. 5c shows a weak negative correlation between bankfull depth and sediment load. Negative correlation is expected, as, the higher the riverine sediment load, the more sediment is deposited on a delta plain, making a steeper delta plain (Fig. 4) and shallower bankfull depth (Fig. 5c).

4 Discussion

Since avulsion is infrequent and often beyond human monitoring timescales, it is difficult to acquire datasets on avulsion frequency from field studies, except for a few well-documented cases (Moodie et al., 2019; Chamberlain et al., 2018; Pierik et al., 2018; Jerolmack, 2009; Li et al., 2022). The six scenarios modelled in this study have merit in providing a large dataset from which to infer the process controls over avulsion timescales.

4.1 Investigating variables controlling avulsion timescales

In this study, to assess the impact of sediment input from the delta's catchment, only the initial alluvial slope (S_alluvial) was varied. All the other measured variables are from these experiments, in which deltas were allowed to self-form. Our experimental setup determines that initial alluvial upstream slope controls how much sediment enters the delta plain (r=0.8; Fig. 4, Table 3). A steeper initial alluvial slope (S_alluvial) has a higher transport capacity and so transports a greater sediment load (Q_s), which is then deposited on the delta plain. Higher sediment load results in a higher vertical aggradation rate (v_a) in the distributary channel (Chadwick et al., 2020), which elevates the distributary channel floor above its surrounding floodplain. The higher the distributary channel floor relative to the floodplain, the easier it is for an avulsion to occur (Jerolmack and Mohrig, 2007; Mohrig et al., 2000). While sediment load is controlled by the steepness of the initial alluvial slope transporting the sediment, we propose that the avulsion timescale in our model is influenced by the amount of sediment deposited in a delta plain. Other investigations that have also found that sediment mass balance is the primary control on avulsion timescales, including a radially averaged model (Muto, 2001; Muto and Steel, 1997), a channel-averaged model (Reitz et al., 2010), and backwater-scaled models (Chadwick et al., 2019; Moodie et al., 2019).

4.2 Comparison with analytical solutions and natural deltas

Figure 6 shows avulsion-timescale scaling relationships observed from our model compared to natural deltas, a physical laboratory fan delta, and analytical solutions. Generally, the pattern and magnitude observed in our numerical model are in good agreement with natural and physical laboratory deltas. By bringing natural and physical laboratory deltas that have a larger range of measured variables (e.g. S_topset, h_c, B, Q_s), clearer correlations between measured variables and avulsion timescale can now be observed. Topset slope negatively correlates with avulsion timescale, showing a steeper delta plain will experience a faster avulsion (Fig. 6a). Positive correlations between bankfull depth, delta lobe width, avulsion length, and avulsion timescale are also consistent between our model and the natural and physical laboratory deltas (Fig. 6b–d). This demonstrates that avulsion takes longer in larger deltas. However, channel width at avulsion, sediment load, and avulsion timescale do not show clear patterns (Fig. 6e, f). Additionally, analytical solutions indicate underestimation of avulsion timescale by several orders of magnitude in comparison to our model and to natural and physical laboratory deltas.

A strong negative correlation between avulsion timescale (T_{a empirical}) and delta topset slope (S_topset) shown in Fig. 6a allows a parallel comparison to be made with avulsion controls in comparable alluvial sedimentary environments (e.g. alluvial fans and fan deltas). Diverse alluvial fan and fan delta experiments have shown that avulsion timescale is influenced by fan-channel gradient (equivalent to delta topset slope in our model) (Schumm et al., 1987; Whipple et al., 1998; Van Dijk et al., 2012; Leenman and Eaton, 2021). Fan-channel slope is dependent on sediment flux (Whipple et al., 1998; Bagnold, 1966) and sediment flux also influences the avulsion timescale as shown earlier in our model. Consequently, avulsion timescale also has a strong correlation with the fan-channel slope or delta topset slope. However, we propose that delta topset slope is a causal effect of the amount of sediment fed into a delta plain in our model. As the initial alluvial slope (S_alluvial) controls the sediment load (Q_s) feeding a delta plain which in turn determines S_topset, we argue that S_alluvial plays a more fundamental role than S_topset in influencing the avulsion timescale (T_{a empirical}) observed in our model.

Table 3 Average sediment load and initial alluvial slope for each scenario.

Figure 5 Ordinary least-squares log–log regressions and power function equations between avulsion timescale and sediment load (a), sediment load and initial alluvial slope (b), and bankfull depth and sediment load (c) measured from the model (N=62) plotted on a log–log scale.

We consider that results from our analytical–numerical model and natural–physical laboratory deltas are directly comparable, but care is needed in their interpretation. The analytical calculations are for fixed values of input variables, and field data are snapshots assumed to represent equilibria. Conditions change during our numerical model runs, and topset slopes reduce through time (Fig. S4a), which would be expected to lead to an increase in avulsion timescale as the delta grows. This results from gentler topset slopes having reduced transport capacities, which reduces the sediment flux, so reducing the in-channel aggradation rate as explained before (Fig. 6a). Similarly, as bankfull depth is also calculated based on the topset slope value (Eq. 1), disagreement between analytical and numerical model results is expected (Fig. 6b). Moreover, delta lobe width (B) and avulsion length (L_a) in the original analytical solution are assumed to be constant with B=40B_c and $L_{\mathrm{a}}=\mathrm{0.2}{L}_{\mathrm{b}}-\mathrm{0.5}{L}_{\mathrm{b}}$ (Chadwick et al., 2020). As we found that delta lobe width (B) and avulsion length (L_a) grow through time (Fig. S4b, c) in our numerical model, analytical model assumptions lead to this disagreement (Fig. 6c, d). Compilation of our numerical model with natural and physical laboratory deltas shows that avulsion becomes less frequent in a larger delta (Fig. 6c, d).

Figure 6 Relationships between avulsion timescales and independent variables: (a) topset slope, (b) bankfull depth, (c) delta lobe width, (d) avulsion length, (e) channel width at avulsion, and (f) sediment load. The plots show model values (T_{a empirical}) as solid black circles. Grey dots and bars are results from the analytical equations using three avulsion threshold H^* values ($T_{\mathrm{a}{H}^{\ast }=\mathrm{1.4}}$, $T_{\mathrm{a}{H}^{\ast }=\mathrm{0.5}}$, $T_{\mathrm{a}{H}^{\ast }=\mathrm{0.2}})$. Diamonds are results from natural and laboratory deltas: grey diamonds indicate examples for which no information about their relative sea-level changes is available; purple diamonds are for deltas with relative sea-level rise (RSLR; mm yr⁻¹) colour-coded as shown. Data from natural deltas and the laboratory experiment are available in Table S2.

The negative correlation between Q_s and T_a from our numerical model (Fig. 5a) deviates from empirical data gathered from natural deltas (Fig. 6f). Natural deltas have relatively larger delta plain sizes than our numerical models represented by avulsion length, bankfull depth, and delta lobe width values (Fig. 6b–d). Assuming no sediment bypass, we normalized sediment discharge to delta size by dividing Q_s by (L_a.B/2) to obtain a vertical aggradation rate, η (m s⁻¹). Vertical aggradation rate is then normalized by dividing it by shear velocity, $u^{\ast }=\sqrt{g{h}_{\mathrm{c}}{S}_{\mathrm{topset}}}$ (m s⁻¹), as we are looking at a vertical aggradation rate related to transport capacity of flow, taken as a measure of the ability to remove sediment from a delta plain with no subsidence, η^* (–). Using this dimensionless aggradation rate, Fig. 7a shows a non-linear negative correlation between T_{a empirical} and η^* from both natural deltas and our models, supporting our argument on sediment discharge influence on avulsion timescale. The best-fit equation for this relationship is $T_{\mathrm{a}\mathrm{empirical}}=-\mathrm{231.5}$ ln(${\mathit{\eta }}^{*})-\mathrm{3573.6}$ (R²=0.71; p<0.05). At ${\mathit{\eta }}^{*}={\mathrm{2.10}}^{-\mathrm{7}}$, the avulsion timescale becomes 0 (Fig. 7a), implying that, at very high aggradation rates with correspondingly high sediment loads, avulsion will effectively occur instantaneously on a delta plain, as the sediment will keep accumulating with very little being removed. Consequently, the flow becomes a random sheet of constant deposition, channel switching, and instantaneous avulsion.

On the other hand, the range of initial alluvial slopes used in the model, although based on observed slopes from 105 global deltas, does not cover the entire range of Q_s inputs to natural deltas (Fig. 6f). Figure 7b shows that a significant increase in Q_s (and a resultant decrease in avulsion timescale) would require the initial alluvial slope to be 6 times the initial downstream slope (Tables 1, 3). Further simulations, using a different initial downstream slope value or more varied slope ratios, could enable a larger range of Q_s to be covered and hypothetically enable a closer agreement between natural deltas and numerical models (Fig. 6f).

Figure 7 (a) Relationship between avulsion timescale (T_a) and dimensionless aggradation rate (η^*) observed from our natural deltas and our models. Regression lines are shown in blue. (b) Boxplots showing distribution of sediment load (Q_s) for each scenario.

The bathymetry of the basin is defined at the beginning of our model runs (Fig. 1). The bathymetry adjusts during a model run to reach an equilibrium profile at a simulation time step > 3–6 d (Video S1). In comparison to other models, in which “reference profiles” and their evolution are defined by making assumptions (e.g. floodplain deposition is assumed to be equal to riverbed aggradation) (Chadwick et al., 2019; Moodie et al., 2019; Moran et al., 2017; Ratliff et al., 2018; Edmonds et al., 2022), our approach directly evaluates avulsion frequency and location that emerge from the physics of self-forming delta lobes. Avulsions in our models consistently arise from channel superelevation (Fig. 1c), consistent with our previous global empirical study (Prasojo et al., 2022).

The ratios between alluvial and topset slope gathered from natural deltas (Table 1) are assumed to be in dynamic equilibrium with the environmental conditions in our study. Even though most modern global deltas have developed since the early Holocene, inevitable natural and anthropogenic changes in boundary conditions, such as changes in sediment or water discharge and local sea-level rise, will have happened, and remain ongoing, during this period (Stanley and Warne, 1994). The slope data derived from remote sensing only represent a single snapshot of a slope that is constantly evolving towards equilibrium with boundary conditions. Consequently, our measured topset slopes from natural deltas may represent a transient condition rather than being equilibrium values. Repeating the work with different initial topset slope presents an opportunity for further investigation.

The avulsion timescales calculated for natural deltas (T_a natural) do not correlate well with the rate of relative sea-level rise (RSLR) (Figs. 6a–f and S5) (Chadwick et al., 2020). As our deltas are self-formed and evolve throughout the simulations, avulsion and backwater lengths (as a function of topset slope; Eq. 1) also grow linearly, and each scenario has a unique avulsion–backwater length ratio (Fig. S6), rather than being constant as previously observed (Ganti et al., 2016b; Chatanantavet et al., 2012). Despite downstream control (e.g. sea-level rise (SLR), backwater length) not being systematically tested in our modelling, empirical correlations between avulsion timescale–SLR (Fig. S5) and avulsion–backwater length (Fig. S6) observed from our numerical models and natural deltas may support our upstream control on delta avulsion hypothesis. We observe that the frequency of avulsions may be unaffected by sea-level rise, as also found in an earlier numerical model study (Ratliff et al., 2021) and a global empirical study (Colombera and Mountney, 2023). Consequently, we propose that avulsion frequency and location are dominated by upstream forcing (i.e. alluvial slope or catchment sediment supply; Figs. 4–7a) (Prasojo et al., 2022) rather than downstream forcing by sea-level rise or backwater effects (Fig. S6) (Chadwick et al., 2020; Chatanantavet et al., 2012; Ganti et al., 2016b). Avulsions in the proximal parts of deltas are controlled by upstream forcing, with RSLR and backwater hydraulics determining avulsion frequency only in distal zones.

Previous literature on the relationship between the frequency of avulsion and sea-level rise is somewhat equivocal. A field study conducted on the Mitchell River delta, Australia, found that avulsion frequency increases with sea-level fall (Lane et al., 2017). Numerical model results suggest that avulsions on the Mississippi (faster) and Trinity (slower) rivers showed different responses to Holocene sea-level rise even though they are geographically adjacent (Chatanantavet et al., 2012; Moran et al., 2017). An example during sea-level fall from the Goose River delta, Canada, shows that avulsion frequency remained constant during this base-level adjustment (Nijhuis et al., 2015). In contrast, avulsion frequency in the Rhine–Meuse delta, the Netherlands, increased during the Holocene sea-level rise period (Törnqvist, 1994), possibly due to aggradation rate (v_a) being controlled by RSLR. Our experiments do not address this issue, and we propose that further investigations combining numerical and flume experiments that are based on observations from natural deltas may aid the resolution of this debate.

4.3 Implications for delta management

Our modelling results advance our understanding of how sediment input from the catchment regulates the timing of avulsions in river deltas. The complex hydraulic and sediment transport processes that lead to the correlation between alluvial slope and avulsion timescale are linked to sediment load, the rate of in-channel aggradation, and how rapidly channels become perched. Consequently, with the increase in anthropogenic forcing both directly within river deltas and throughout upstream catchment areas (Best, 2019; Darby et al., 2015; Dunn et al., 2019; Hackney et al., 2020), river delta managers can use sediment load management to reduce the risk of avulsion-driven flooding.

However, finding a perfect balance between reducing avulsion frequency, maintaining the sediment load required to nourish delta environments and to hinder deltas' risk from subsidence and coastal erosion, is challenging. In some locations, deforestation that increases sediment supply is responsible for ∼ 25 % of net land gain on global deltas, which may also hasten future avulsions (Nienhuis et al., 2020). Conversely, river impoundment is responsible for a more than 50 % reduction in sediment delivery to the global ocean since 1950, collectively leading to a loss in delta land of 127±8.3 km² annually over the 30 years from 1980–2010 (Nienhuis et al., 2020). This declining sediment input not only poses threats to the long-term sustainability of deltas but also renders them susceptible to adverse effects from rising sea levels and ecological degradation due to sediment starvation and saltwater ingress (Jordan et al., 2019). Therefore, gaining insights into the distribution patterns and quantities of sediments in deltas is imperative to ensure their continued sustainability.

4.4 Next steps

An important extension of this modelling work is to have more varied S_topset : S_alluvial ratios, water discharge (Q), and sediment load (Q_s), as variability in these may affect the geomorphic processes controlling avulsion timescale. Moreover, investigation using scenarios based on sets of geological variables that control slopes, rather than varying the slopes themselves, will be regarded as valuable. Multi-temporal observation of well-studied natural river deltas, such as the Yellow (Moodie et al., 2019), Mississippi (Chamberlain et al., 2018), or Rhine–Meuse (Pierik et al., 2018) deltas, could then be used to validate model results. On the other hand, incorporating other variables, such as grain size and sediment cohesion, forcing through sea-level rise and subsidence, varying basin geometry, and adding vegetation that controls crevassing and consequently increases avulsion timescale in future numerical modelling, should be considered and are already in other studies (Nienhuis et al., 2018a; Pierik et al., 2023; Sanks et al., 2022; Wolinsky et al., 2010). In particular, considering the importance of projected global sea-level changes and the variability in results reported in the literature, a better understanding of sea-level rise impacts on delta avulsion is needed.

We have used a simplified modelling approach and have isolated one predictor variable while holding other factors constant. Observations of the processes and evolution in the numerical deltas show the complexity of hydraulic and morphodynamic processes across delta plains. Future work will need to address this complexity by answering the following questions. (a) How does the forcing studied here (alluvial slope and consequent sediment input) interact with a combination of other factors (e.g. sea level, wave and tidal regimes, and anthropogenic effects)? (b) How do the other controls (e.g. Q_s, Q, riverbank material, vegetation) in river deltas influence avulsion timescales? (c) How might these avulsion signals be preserved or shredded in the rock record?

5 Conclusions

We conducted a suite of numerical morphodynamic modelling experiments with variable river alluvial slopes (from $\mathrm{1.13}\times {\mathrm{10}}^{-\mathrm{4}}$ to $\mathrm{3.04}\times {\mathrm{10}}^{-\mathrm{3}}$) to understand the controls over avulsion location and timescale in a river delta. Sediment load, which in this study was controlled by the imposed alluvial slope upstream of a delta plain, controls the avulsion timescale. Mechanistically, when channel width is constant, a steeper alluvial slope has a greater sediment transport capacity. Consequently, a steeper alluvial slope may transport a larger sediment load to be delivered to a delta plain. This induces accelerated vertical aggradation that causes more frequent avulsion in our model. Integrating our numerical modelling results with analytical, natural, and physical laboratory deltas supports the hypothesis of upstream forcing influencing delta avulsion timescale and location, rather than downstream influence by backwater length or sea-level rise. However, our model has several limitations, such as dynamic equilibrium assumption from our alluvial–topset slope ratios, homogeneous initial topset slope adopted from Mississippi Delta, and constant discharge applied in the model, along with factors that were ignored (e.g. tectonics). Alluvial slope exerts a strong control on the timing of avulsion in river deltas and thus directly impacts (i) coastal and inland hazards on highly populated river deltas and (ii) the interpretation of deltaic sedimentary deposits.