Marsh-induced backwater: the influence of non-fluvial sedimentation on a delta's channel morphology and kinematics

We investigate the interaction of fluvial and non-fluvial sedimentation on the channel morphology and kinematics of an experimental river delta. We compare two deltas: one that evolved with a proxy for non-fluvial (“marsh”) sedimentation (treatment experiment) and one that evolved without the proxy (control). We show that the addition of the non-fluvial sediment proxy alters the delta's channel morphology and kinematics. Notably, the flow outside the channels is significantly reduced in the treatment experiment, and the channels are deeper (as a function of radial distance from the source) and longer. We also find that both the control and treatment channels narrow as they approach the shoreline, though the narrowing is more pronounced in the control compared to the treatment. Interestingly, the channel beds in the treatment experiment often exist below sea level in the terrestrial portion of the delta top, creating a ∼ 0.7 m reach of steady, non-uniform backwater flow. However, in the control experiment, the channel beds generally exist at or above relative sea level, creating channel movement resembling morphodynamic backwater kinematics and topographic flow expansions. Differences between channel and far-field aggradation produce a longer channel in-filling timescale for the treatment compared to the control, suggesting that the channel avulsions triggered by a peak in channel sedimentation occur less frequently in the treatment experiment. Despite this difference, the basin-wide timescale of lateral channel mobility remains similar. Ultimately, non-fluvial sedimentation on the delta top plays a key role in the channel morphology and kinematics of an experimental river delta, producing channels which are more analogous to channels in global river deltas and which cannot be produced solely by increasing cohesion in an experimental river delta.

1 Introduction

River deltas are dynamic systems that accumulate sediment through the interaction of many processes occurring in the distributary channels, floodplain, and basin. The formation and movement of channels (channel kinematics) distribute fluvial sediment across the delta plain, leading to both active depocenters and quiescent regions that change location through time (i.e., river avulsions; e.g., Coleman, 1988; Straub et al., 2009; Hoyal and Sheets, 2009; Chadwick et al., 2020). The active depocenters are dominated by fluvial sedimentation, which includes sediment deposition in mouth bars, channels, and the delta front, as well as finer-grained overbank deposition adjacent to the river occurring during floods (e.g., Coleman, 1988; Paola et al., 2011; Khan et al., 2013). The passive areas aggrade through fine-grained mineral sediment delivered from the ocean through waves, tides, and storms (e.g., Smith et al., 2015; Sanks et al., 2020), primary production of organic material via wetland vegetation (e.g., Nyman et al., 2006; Holmquist et al., 2018, 2021; Kelsall et al., 2023), and/or the redistribution of sediment from marsh edge erosion (e.g., Hopkinson et al., 2018; Valentine et al., 2023), the combination of which we refer to as non-fluvial (“marsh”) sedimentation.

Understanding and managing sediment accretion in rivers, floodplains, and adjacent wetlands remains a key strategy to mitigate the impacts of relative sea level rise (RSLR) worldwide. For example, Louisiana has planned and implemented both engineered sediment diversions (e.g., the recently approved mid-Barataria sediment diversion) and artificial marsh creation projects (e.g., CPRA, 2017; Nittrouer et al., 2012; White et al., 2019) to reconnect the river to its wetlands and restore the coastal land area (e.g., Elsey-Quirk et al., 2019; Peyronnin et al., 2017; Xu et al., 2019). Another strategy is engineered river avulsions, in which the avulsion location can be optimized if the location of the maximum channel sedimentation is known (Moodie and Nittrouer, 2021). On the short-term, the success of these mitigation strategies depends on the interaction of channel and wetland aggradation (e.g., Kirwan and Megonigal, 2013; Paola et al., 2011), which can be maximized under certain conditions (Esposito et al., 2017). This is because, in order for these ecosystems to succeed on short (management) timescales, they need to accumulate both organic sediment produced in situ and mineral sediment from channels or offshore processes (e.g., Sanks et al., 2020; Elsey-Quirk et al., 2019; Kelsall et al., 2023). However, the implications of wetland–channel coupling over timescales longer than avulsions (10–10³ years, depending on the river delta) remain unclear. As such, understanding the equilibrium morphology of a coupled wetland–delta system can enhance short- and long-term management decisions.

The deposition of fluvial sediment is controlled by the proximity to channels in a river delta, and the life of a delta depends on the intrinsic ability for a river to avulse (Slingerland and Smith, 2004). Fluvial sedimentation typically occurs via topographic flow expansions in physical delta experiments (e.g., Hoyal and Sheets, 2009; Chatanantavet and Lamb, 2014; Sittoni et al., 2014; Shaw et al., 2018). As the river enters the basin, the channels lose confinement, which triggers rapid deposition at the mouth bar (e.g., Edmonds and Slingerland, 2007; Sittoni et al., 2014; Törnqvist and Bridge, 2002). The rapid deposition leads to an upstream migration of the depocenter, which is referred to as the “morphodynamic backwater effect” (Hoyal and Sheets, 2009; i.e., “backfilling” in the alluvial fan literature; Bull, 1968). Eventually, this process leads to a channel avulsion, which is thought to occur at some distance upstream of the shoreline and is related to channel superelevation (e.g., Edmonds et al., 2009; Shaw et al., 2021; Mohrig et al., 2000; Hajek and Wolinsky, 2012; Ganti et al., 2016b; Jobe et al., 2020). The location of avulsions appears to occur preferentially at the backwater length scale (Jerolmack and Swenson, 2007), which is defined as the channel depth divided by the energy slope for low Froude number systems (e.g., Shaw and McElroy, 2016). This length scale serves as an estimate of the distance between a channel mouth and the location where the channel bed of a river drops below sea level (Paola and Mohrig, 1996). Morphodynamic models of the channel deposition suggest that avulsion locations depend on variable discharge of a river to move the location of maximum aggradation upstream through scour during floods (e.g., Nittrouer et al., 2011, 2012; Chatanantavet et al., 2012; Ganti et al., 2016b; Chadwick et al., 2019). It is argued that when fluctuating flows are not present, the location of the maximum deposition is inevitably at the channel mouth, producing topographic flow expansions (or morphodynamic backwater effects) that typically dominate the channel kinematics of physical delta experiments.

Avulsion locations scale with hydrodynamic backwater length, but there appear to be multiple possible causes for this scaling. The backwater length, or the length of the river that is influenced by offshore processes, is thought to be controlled by complex hydrodynamics of rivers (e.g., Lamb et al., 2012; Chatanantavet et al., 2012; Nittrouer et al., 2011). More recently, Ratliff et al. (2021) and Prasojo et al. (2022) have shown that avulsion locations on a river delta can be explained without backwater hydrodynamics and are a function of a break in slope in the channel (i.e., morphology of the channel). In other words, even though the avulsion location scales with the hydrodynamic backwater length, Ratliff et al. (2021) suggest that avulsion locations and subsequent channel movement can be explained simply through a change in geometry without the need for complex backwater hydrodynamics. Similar to Ratliff et al. (2021), Moodie et al. (2019) showed that a decrease in channel slope (i.e., a change in channel geometry) triggers lobe progradation, subsequently moving the location of maximum channel aggradation (and the avulsion location) further upstream. Relatedly, Sanks et al. (2022) showed that non-fluvial aggradation near the coastline can significantly rearrange the depositional patterns on experimental deltas, thus influencing delta hypsometry and floodplain geometry. As such, non-fluvial deposition may control channel kinematics as well.

Despite these different mechanisms for avulsions, the impact of non-fluvial sedimentation on the channel kinematics of deltaic systems remains largely unknown (Paola et al., 2011). While channel morphology and kinematics of river deltas in experimental (e.g., Hoyal and Sheets, 2009; Li et al., 2017; Straub et al., 2013; Barefoot et al., 2021; Carlson et al., 2018), numerical (e.g., Edmonds et al., 2009; Caldwell and Edmonds, 2014; Lauzon and Murray, 2018), and field settings (e.g., Shaw et al., 2016; Wilson and Goodbred, 2015; Aslan et al., 2005; Carlson et al., 2021) are well documented, wetland dynamics are rarely seen as being coupled to, or a driver for, channel dynamics. In numerical and experimental settings, added sediment cohesion is often used to represent the influence of vegetation on the delta top for its ability to increase critical shear stresses and reduce erosion (e.g., Hoyal and Sheets, 2009; Edmonds and Slingerland, 2010; Li et al., 2017). Another numerical study shows that there is an optimal vegetation height for inorganic sediment trapping on deltaic marsh platforms (Nardin and Edmonds, 2014). Similarly, Piliouras et al. (2017) used alfalfa to simulate vegetation in an experimental delta, which showed enhanced sediment trapping and increased delta-top roughness. Because the alfalfa increased the delta-top slope, they found that the results of the experiments were only applicable to fan deltas and not the low-sloping coastal deltas of interest here. While valuable, these previous studies leave out a critical property of wetlands, which is the ability to accumulate non-fluvially delivered sediment in situ via the primary production of vegetation (i.e., accumulation of above- and belowground biomass; e.g., Kosters et al., 1987; Morris et al., 2002; Mudd et al., 2009) and the accumulation of fine-grained mineral sediment (muds), presumably delivered from the ocean (e.g., Leonardi et al., 2021; Sanks et al., 2020). In contrast to this body of work, Sanks et al. (2022) and Zapp et al. (2022) show that non-fluvial sedimentation can act independently of cohesion by changing the distribution of elevations, delta-top slope, sediment partitioning, and shallow compaction rates.

Figure 1 (a) Aerial image from the treatment experiment adapted from Sanks et al. (2022) showing the entrance channel, which delivers water and fluvial sediment to the delta top, and the marsh sediment dispenser, which delivers non-fluvial sediment to the delta top. The brown sediment is the kaolinite marsh proxy, which forms a platform near sea level. (b) Conceptual diagram of a delta–wetland system adapted from NASA (2014). In field systems, the platform aggrades through a combination of fluvially delivered fine-grained sediment, muds and organics delivered from offshore, and organic sediment produced in situ by vegetation. In the treatment experiment, the platform encompasses the area of the delta between 5 and −9 mm relative to sea level and aggrades through our simplified elevation-based model. The fluvial deposit exists above 5 mm relative to sea level in the treatment experiment and aggrades solely through the deposition of fluvially delivered sediment.

To test the control of non-fluvial sedimentation on channel kinematics, we separate the broad variation in delta-top facies into “non-fluvial” (e.g., wetlands, interdistributary bays, and tidal flats) and “fluvial” (e.g., upper delta plain, levees, mouth bars, and prodelta) components (Bhattacharya, 2006). We define non-fluvial sedimentation as the spatially extensive, persistent, fine-grained, and compactible deposition that occurs in the quiescent, protected regions of deltas most commonly associated with wetland platforms (Fig. 1). In contrast, the fluvial sediment that is the default on experimental deltas is coarser grained, less compactible, and not persistent. If a small amount of non-fluvial sedimentation can significantly rearrange the delta mass balance (Sanks et al., 2022), then it is possible for it to affect channel dynamics as well. We investigate this hypothesis here.

2 Methods

2.1 Experimental setup and data

Despite the differences in the scale of field and experimental systems, physical experiments approximate conditions and processes that occur in global systems (Paola et al., 2009). Thus, physical delta experiments are an ideal way to study the influence of non-fluvial sedimentation on channel morphology and kinematics. We compare two experiments run under the same boundary conditions (Table 1), except that one experiment has a proxy for non-fluvial (marsh) sedimentation (treatment) and one does not (control). The addition of the kaolinite marsh proxy accounts for ∼ 15  % of the final delta volume and ∼ 8  % of the mass in the treatment experiment (Sanks et al., 2022). Thus, we assume that any statistically significant changes in channel kinematics can be attributed directly to this non-fluvial sedimentation. Both experiments were run for 560 h or ∼ 20 times the compensation timescale to capture autogenics (i.e., self-generated dynamics like channel avulsions) and account for the inherent stochasticity of the systems (Wang et al., 2011). Ultimately, the systems accumulate about 20 channel depths of stratigraphy, as the channels are ∼ 7 mm deep, and the deltas evolve with a background relative sea level rise rate (RSLR_b) of 140 mm over 560 h. Because the experiments are run in an equilibrium state, we do not focus on temporal dynamics herein, though this is an interesting area for future research.

Table 1 The experimental conditions for both the control (no marsh) and treatment (marsh) experiments used for comparison in this study taken from Sanks et al. (2022).

Elevation relative to sea level (r.s.l.) is a primary control on the deposition of mud in tidal flats and wetland platforms (e.g., Fagherazzi et al., 2006) and organic sediment production in wetlands (e.g., Morris et al., 2002). Thus, an elevation-based model is representative of all non-fluvial sedimentation that occurs in deltaic coastal environments. Note that we neglect the drag of vegetation in this study, as we are interested solely in the impact of additional mass that accumulates in global deltas on the long-term dynamics of the system. Although wetland sedimentation can be a complex process, we adapt a model tying primary production in salt marshes to elevation relative to mean high tide from Morris et al. (2002) and simplify this model to produce three distinct regions of non-fluvial deposition, namely unstable (−9 to −5 mm r.s.l.), maximum production (−5 to 0 mm r.s.l.), and stable (0 to 5 mm r.s.l.; Fig. 1b). Because the experiments evolve with the absence of tidal processes, we scaled the marsh window to the relative sea level instead of mean high tide. The marsh window spans 14 mm in elevation, which is equivalent to about one experimental channel depth. The marsh window elevation band encompassed 30 % of the planform area in the control experiment. The non-fluvial sediment is Edgar Plastic Kaolin (a type of kaolinite clay with a mean grain size of 1.36 µm) that we deposit from above using a sieve (see the Methods section in Sanks et al., 2022; Fig. 1a). The river sediment mixture has a range of grain sizes from 1 to 1000 µm, with a mean of 67 µm, and contains a polymer for added cohesion (Straub et al., 2015). Deposition rates are determined using a hexagonal grid (7.5 cm sides) to find the average elevation inside each hexagon. If the average elevation falls within one of the three regions, we deposit either 3.4 g (maximum production; accumulates 1RSLR_b) or 1.7 g (stable and unstable; accumulates 0.5RSLR_b) of kaolinite every 2 h in that hexagon. The deposition rate was calibrated based on the assumption that kaolinite deposited in water has a porosity of 90 %. Non-fluvial deposition occurred while the experiments were running, so that any kaolinite deposited in the channels would be immediately washed out and transported offshore.

To analyze the channel properties and kinematics from both experiments, we use dry and wet lidar scans, binary channel maps, and binary flow maps (Appendix A). We collect dry lidar scans every hour for the control and every 2 h for the treatment while the experiments are paused. Because non-fluvial deposition occurs every 2 h in the treatment experiment, this controls the frequency of dry lidar scans. As explained later, all sedimentation rate analyses are conducted on a 2 h timescale for both experiments to avoid the Saddler effect. The dry lidar scans have 5×5 mm horizontal resolution and a sub-millimeter vertical resolution and are used solely for elevation-based analyses. We collect wet lidar scans every hour for both experiments while they are running. The wet lidar scans are used for flow mapping and also have 5×5 mm horizontal resolution. We create binary channel maps by hand mapping the channels using hourly aerial imagery (generated from the wet scans; Fig. A1) because simple color thresholds are unable to differentiate between channelized and non-channelized flow. However, a simple color threshold can capture total flow on the delta top (Fig. 3c, d), so we use this method to create total flow maps. As such, overbank flow maps can be produced by differencing the total and channel flow maps. While overbank maps are sensitive to thresholding, the final results are insensitive to threshold choice. We note that due to artifacts from the moving cart during marsh deposition, channels were obstructed in some aerial images (n=73) in the treatment experiment.

2.2 Channel properties

We analyze the treatment and control experiments for differences in channel properties (Table 2), including but not limited to channel depth (H_c), trunk channel depth (H_tc), channel planform area (A_c), channel length (L_c), and trunk channel width (W_tc). Channel properties are calculated for each channel segment within a radial transect. Channel depths are calculated every 5 mm (width of 1 lidar pixel) from 0 mm (the apex) to 3100 mm (max channel length) from the apex, while all other channel properties are calculated every 50 mm from 0 mm (the apex) to 3100 mm from the apex because those mean properties are not sensitive to bin width. All trunk channel depths are calculated by taking the deepest channel in each transect, and all trunk channel widths are calculated by taking the widest channel in each transect. In any context that we do not specify that we are analyzing only the trunk channel, we are analyzing all channels on the delta top (trunk and non-trunk channels).

To calculate channel depths, we use a square channel buffer of 20 mm (4 pixels) on all sides of the channel maps to ensure that the levee crests are captured in the measurement. H_c (mm) is then given by

$$\begin{array}{}\text{(1)}& H_{\mathrm{c}}=max\left(z_{\mathrm{c}}\right)-min\left(z_{\mathrm{c}}\right),\end{array}$$

where max(z_c) is the maximum channel elevation for each channel within the radial transect (i.e., levee crest; Fig. 2c, d), and min(z_c) is the minimum channel elevation for each corresponding channel in the radial transect (i.e., channel thalweg; Fig. 2c, d). We use the levee crest (maximum channel elevation) to compute the channel depth because the flow is almost always at the levee crest. Because each transect often has multiple channel segments and thus multiple channel depths, we compute a trunk channel depth (H_tc; mm) for each radial segment through time by

$$\begin{array}{}\text{(2)}& H_{\mathrm{tc}}=max\left(H_{\mathrm{c}}\right).\end{array}$$

Now, every 5 mm from the apex, we have one channel depth (i.e., the trunk channel depth) for each time step. We use these data to compute the mean channel depth (through time) for H_tc as a function of distance from the apex. We also determine the basin-wide $\stackrel{\mathrm{‾}}{H_{\mathrm{c}}}$, which is simply the mean of all (trunk and non-trunk) channel depths through space and time.

The basin-wide channel area (A_cb; cm²) is the time-averaged channel area (A_c). A_c is simply the sum of all channel pixels on the delta top multiplied by 0.25 cm² (the area of 1 pixel). Mean radial channel fraction $\stackrel{\mathrm{‾}}{f_{\mathrm{Ac}}}$ is then given by

$$\begin{array}{}\text{(3)}& \stackrel{\mathrm{‾}}{f_{\mathrm{Ac}}}={\displaystyle \frac{A_{\mathrm{cr}}}{A_{\mathrm{r}}}},\end{array}$$

where A_cr is the time-averaged channel area of the radial transect, and A_r is the area of the radial transect.

The channel length (L_c) is assumed to be the straight-line distance from the entrance channel to the most distal channel pixel at each time step. We note that this is an underestimate for the length of any sinuous channels in the experiments, but a visual inspection shows nearly straight channels, with no systematic variation between experiments. W_tc is the trunk channel width (cm) and is given by

$$\begin{array}{}\text{(4)}& W_{\mathrm{tc}}=max(A_{\mathrm{cr}}/L_{\mathrm{r}}),\end{array}$$

where A_cr is the area of each channel in the 5 cm radial bin of interest (cm²), and L_r is 5 cm, which is the length of each radial transect. We note that this is likely an overestimate for the last radial transect (i.e., the channel tip) because the channelized portion in the last radial transect likely does not encompass the entire length of the radial transect.

Figure 2 (a) Lidar map of the control delta at run hour 81, showing elevations relative to sea level. The black line represents the transect A-A^′, which is 1 m from the apex. See Fig. 8a for an aerial image. (b) Lidar map of the treatment delta at run hour 221, showing elevations relative to sea level. The black line represents the transect B-B^′, which is 1 m from the apex. See Fig. 8b for an aerial image. (c) Elevation relative to sea level versus the distance along the radial transect A-A^′. (d) Elevation relative to sea level versus the distance along the radial transect B-B^′.

We are also interested in the hydrodynamic backwater effects present in both experiments. We define the hydrodynamic backwater length as the longest continuous length of the channel reach, where the channel bed is at or below sea level (Chatanantavet et al., 2012). To calculate the length of the backwater reach, we use dry lidar scans and channel maps at 1 h (control) and 2 h (treatment) intervals (n=560 for the control and n=265 for the treatment, due to artifacts rendering 15 time steps unusable). We determine the location of the beginning of the longest continuous backwater reach for each time step, and this location is referred to as the backwater point (Fig. 2b). If the channel tip is above sea level or there is no radial channel transect with at least 16 % of channel bed elevations at or below sea level, then there is effectively no backwater reach, and we assume that the backwater length is zero. Although the absolute backwater length is somewhat sensitive to this threshold choice, the presence of a significant backwater length in the treatment experiment as compared to the control is not.

2.3 Channel kinematics

Various metrics were used to compare channel kinematics between the experiments. The delta-top area used to quantify lateral channel movement is defined as the area above sea level for at least 50 % of the experiment (p_land; Fig. B1), consistent with Li et al. (2017). We compute the fraction of the delta (delta area is p_land; m²) that has not been visited by a channel through time, which can be used to estimate the lateral channel mobility of the system (Li et al., 2017). This metric (f_mob; unitless) is defined as

$$\begin{array}{}\text{(5)}& f_{\mathrm{mob}}={\displaystyle \frac{\sum _i^n{A}_{\mathrm{c}}}{p_{\mathrm{land}}}},\end{array}$$

where i is the initial hour, and n is every time step from 1 to 560 (e.g., if using a time step n>i, then the numerator is the cumulative channel area from hours i through n; A_c; m²). The unvisited area of the experiments decays exponentially through time. To describe this, we compute the e-folding lateral mobility timescale (T_mob; h). T_mob is the inverse of the exponent from the function that describes the lateral channel mobility decay of the system (Li et al., 2017). Radial lateral mobility is also calculated every 50 mm from 0 mm (the apex) to 3100 mm from the apex. In this case, A_c is the channelized area in the radial transect, and p_land is the area of the radial transect. Similarly, we calculate the planform overlap decorrelation metric described in Wickert et al. (2013) as an independent check on the lateral mobility of the system. This metric describes the time it takes for subsequent channel maps to decorrelate. The planform overlap is described in Appendix B (see also Fig. B5). Both the lateral channel mobility and planform overlap indirectly describe the avulsion timescale, as they determine how quickly channels move across the delta top. While channels move in other ways besides avulsions, the largest changes occur due to channel avulsions, and thus, the time it takes for channels to visit most of the delta top (lateral mobility) and how long it takes channel maps to decorrelate from each other (planform overlap) are channel movement statistics.

We then calculate channel mobility metrics related to sedimentation patterns on the delta, as channel sedimentation can trigger an avulsion (e.g., Jobe et al., 2020; Ganti et al., 2016b; Moodie and Nittrouer, 2021). The first is the fraction of the delta that is unmodified (f_um; unitless), which is defined as the fraction of the delta that has accumulated at least 1 mm of sediment in the terrestrial delta (p_land) for various periods of time. Similar methods to the lateral mobility are used to determine f_um (Li et al., 2017). The modification timescale (hours) is the inverse of the exponent from the function that describes how long it takes to accumulate 1 mm of sediment in the terrestrial delta.

Next, we calculate the compensation timescale (T_c; h). Compensation is a metric used to describe surface processes in channelized systems and compares long-term sedimentation patterns to general accommodation of the system. In general, compensation describes a channel's inherent tendency to fill low-lying areas on the delta (Straub et al., 2009). T_c refers to how long it takes to accumulate roughly one channel depth of sediment everywhere on the delta top. Again, to capture levee crests in the channel depths, we use a square channel buffer of 20 mm (4 pixels) on all sides of the channel maps. The basin-wide T_c is given by

$$\begin{array}{}\text{(6)}& T_{\mathrm{c}}={\displaystyle \frac{H_{\mathrm{c}}}{V_{\mathrm{t}}}},\end{array}$$

where H_c is the basin-wide mean channel depth (mm), and V_t is 0.25 mm h⁻¹, which is the mean aggradation rate for both experiments. T_c is also calculated radially from the apex using 5 mm bins, where H_c is the mean trunk channel depth (mm) of the radial bin, and V_t is the mean total aggradation rate (mm h⁻¹) of the radial bin. To directly compare aggradation rates between the experiments at the same timescale (Sadler, 1981), we use 2 h lidar difference maps for both experiments.

Last, we determine the channel in-filling timescale (T_f; h). Avulsion locations can be related to peaks in channel aggradation (e.g., Chatanantavet et al., 2012; Chadwick et al., 2019; Moodie and Nittrouer, 2021), so we calculate channel and far-field aggradation rates (mm h⁻¹), allowing us to compare channelized (channel) and non-channelized (far-field) sediment depositional patterns. We again use the buffered channel maps and compute aggradation rates every 2 h. The channel in-filling timescale (T_f; h) is given by

$$\begin{array}{}\text{(7)}& T_{\mathrm{f}}={\displaystyle \frac{H_{\mathrm{c}}}{V_{\mathrm{c}}-V_{\mathrm{ff}}}},\end{array}$$

where H_c is the basin-wide mean channel depth (mm), V_c is the basin-wide channel aggradation rate (mm h⁻¹), and V_ff is the basin-wide far-field aggradation rate (mm h⁻¹). We also compute a channel in-filling timescale radially from the apex every 5 mm. In this case, H_c is the time-averaged trunk channel depth for each radial transect, V_c is the time-averaged channel aggradation rate (mm h⁻¹) for each radial transect, and V_ff is the time-averaged far-field aggradation rate (mm h⁻¹) for each radial transect.

3 Results

3.1 Channel properties

Does the presence of non-fluvial sedimentation influence channel morphology in experimental deltas? A simple comparison reveals considerable differences. Importantly, there is a reduction in the mean fraction of the delta top covered in overbank flow by about one-third in the treatment experiment (0.183 ± 0.122; all error bounds represent 1 standard deviation from the mean, unless otherwise noted) compared to the control (0.489 ± 0.119). Though this decrease in overbank flow does not lead to an increase in planform channel area, it does suggest that the flow is concentrated in channels in the treatment experiment. A loss of the downstream channel confinement is observed, which is associated with an increase in overbank flow in both the control and treatment experiments near the mean shoreline (Fig. 3a and b).

Figure 3 (a) Mean channelized and (b) overbank flow fraction as a function of distance from the apex for the control (blue) and treatment (green), with 1σ standard deviation shown in the shaded polygons. (c) The channelized and overbank flow on the terrestrial delta top for the control experiment at hour 181, which approximates the mean flow conditions during the experiment. (d) The channelized and overbank flow on the terrestrial delta top for the treatment experiment at hour 360, which approximates the mean flow conditions during the experiment.

The basin-wide channel morphology of the control and treatment experiments is different (Table 2). The channels are on average narrower (Fig. 4b) and longer (Fig. 4c) in the treatment experiment compared to the control. Furthermore, the treatment experiment has a greater number of distributary channels (computed automatically through image segmentation) on the delta top (Fig. 4d).

Table 2 Basin-wide channel properties for the control and treatment experiments.

The channels in the treatment and control experiments extend past the mean shoreline (Fig. 4a and c). The mean shoreline in the control experiment is 0.942 ± 0.156 m (blue open circles in all radial figures; the interquartile range (IQR) is 0.874 m [0.486,1.36]), and the mean channel length is 1.51 ± 0.211 m. On average, the channels extend 0.571 ± 0.265 m past the mean shoreline. The mean shoreline in the treatment experiment is 1.11 ± 0.156 m (green open circles in all radial figures; the IQR is 0.889 m [0.651,1.14]), and the mean channel length is 1.90 ± 0.420 m. On average, the channels extend 0.793 ± 0.413 m past the mean shoreline. While the slope break occurs around the median shoreline (∼ 1.5 m from the apex) for the treatment experiment (i.e., the channels extending past the local shoreline (overstepped channels) have a slope break), there is no such break in the control (Fig. 4a), even though the channels extend past the mean shoreline as well. The channel beds in the treatment are often below mean sea level, suggesting the presence of hydrodynamic backwater effects (Fig. 4a). However, this is not the case in the control, as the channels are perched above sea level, suggesting that the channels do not always extend past the local shoreline. The channels in the control experiment are wider near the apex and get narrower as they approach the shoreline. The distal narrowing of the channels is also observed in the treatment experiment, but the variation is channel width as a function of distance from the apex is not as large (Fig. 4b). The combination of more and longer channels in the treatment experiment produces channels with similar planform area as the control (Table 1). The number of channels decreases with radial distance from the apex in both experiments, supporting loss of channel definition (i.e., a loss of self-formed levees) near the shoreline.

Figure 4 Channel morphology for the control (blue) and treatment (green) experiments. (a) The mean channel bed elevation relative to sea level (mm) as a function of radial distance from the apex (m). (b) The mean trunk channel width (cm) as a function of radial distance from the apex (m). (c) Violin plots of the channel length (m) for the control and treatment experiments. (d) The mean number of channels as a function of distance from the apex (m). The shaded blue (control) and green (treatment) areas signify 1σ standard deviation from the mean, and open circles are the mean shoreline positions.

How does the addition of marsh deposition affect the presence of a backwater reach? We show that the treatment experiment has a significant backwater reach that is not observed in the control experiment (Fig. 5). Out of the 560 h of the control experiment, there is no radial transect with at least 16 % of the elevations at or below sea level for 49 % of the experiment (274 h). Of the 286 hours in the control experiment with an observed backwater reach, only 1 % of the channels have backwater lengths >  0.5 m, with a maximum backwater length of 0.57 m (Fig. 5). However, a backwater length is always observed in the treatment experiment for the 265 h where there is a viable channel map. The treatment experiment has a backwater length >  0.5 m 69 % of the time, with a maximum backwater length of 1.91 m (Fig. 5). The average backwater length in the control experiment is 0.0474 m, and it is 0.683 m in the treatment (Table 2; Fig. 5). Both the control and treatment experiments show variability in backwater lengths through time, with peaks roughly every 50 h (∼T_c).

Figure 5 (a) The backwater length (m), measured as the maximum consecutive length of the channel domain that is at or beneath sea level, for the control (blue) and treatment (green) experiments through time. The black crosses represent the times with no channel maps for the treatment experiment. (b) Violin plots of the backwater length for the control (blue) and treatment (green) experiments.

3.2 Channel kinematics

Despite the differences in channel morphology, we observe only subtle differences in the statistics characterizing channel kinematics or the motion of the channels. We first characterize the channel and far-field aggradation rates for both experiments and observe differences (Fig. 6a). The mean channel aggradation rate in the control experiment increases with radial distance from the apex with a large peak near the downstream end, which is characteristic of topographic flow expansions (Fig. 6a). In contrast, the treatment experiment shows a significantly reduced peak in the mean channel aggradation rates near the shoreline (Fig. 6a). The long-term aggradation rate near the mean shoreline for both experiments is 0.25 mm h⁻¹ (Sanks et al., 2022), but the highly mobile channels mean that, at short timescales, channel aggradation rates near the shoreline are more rapid than long-term aggradation rates. The reduced peak in the treatment channels compared to the control channels corresponds to the start of the hydrodynamic backwater zone (Fig. 4a) and the larger number of channels in that zone (Fig. 4b).

Figure 6 (a) The mean channel (solid line) and far-field (dashed line) aggradation rates (mm h⁻¹) in the control (blue) and treatment (green) experiments as a function of radial distance from the apex (m). (b) The mean trunk channel depth (mm) as a function of radial distance from the apex (m). (c) The aggradation difference (mm h⁻¹; i.e., the channel in-filling rate) between far-field and channel aggradation rates as a function of radial distance downstream (m). (d) The channel in-filling timescale (h) as a function of radial distance downstream (m) for locations ≥ 0.1 m from the apex. Please refer to Fig. C1 for the same figure which shows the data that have been cut off from panels (a) and (c) here.

Sediment accumulation rates show that the mean basin-wide channel aggradation rate for the control experiment is 0.425 ± 1.23 and 0.450 ± 1.15 mm h⁻¹ for the treatment. The radially and time-averaged channel aggradation rate for the control experiment is 1.06 ± 0.869 and 0.679 ± 0.190 mm h⁻¹ for the treatment (Fig. 6a). The basin-wide mean far-field aggradation rate for the control is 0.050 ± 2.21 and 0.100 ± 1.38 mm h⁻¹ for the treatment. The radially and time-averaged far-field aggradation rate for the control experiment is 0.277 ± 0.180 and 0.250 ± 0.146 mm h⁻¹ for the treatment (Fig. 6a). Because the basin-wide channel depths are similar for both experiments, this results in a basin-wide channel in-filling timescale (T_f) of 8.70 h for the control and 13.8 h for the treatment experiment (using Eq. 7). It takes about 1.5 times as long for channels to fill in relative to the floodplain in the treatment experiment than the control, suggesting the presence of channel buttressing via the non-fluvial sedimentation in the treatment experiment. Despite similar basin-wide channel depths, the trunk channels in the treatment experiment are deeper than in the control (Fig. 6b), again producing a longer T_f as a function of radial distance from the apex (Fig. 6d).

Both control and treatment experiments exhibited mobile channels with some clear episodes of avulsion and lateral migration (Appendix B). However, the steady presence of shallow unchannelized flow in both experiments made it impossible to map the discrete avulsion locations. Instead, both deltas appeared to maintain overbank flows that kept the “finding phase” (Reitz et al., 2015) of the avulsion cycle active (Figs. B3, B4). No clear distinction in avulsion location or timescale could be established (e.g., Figs. 6, 7, and B5–B8), despite the significant difference in backwater length and aggradation patterns.

The basin-wide lateral mobility is roughly the same in both experiments (Fig. 7a). The basin-wide lateral mobility timescale (T_mob) for the control experiment is 63 and 60 h for the treatment experiment (Fig. 7a). However, we see subtle differences in T_mob as a function of radial distance from the apex (Fig. 7c, d). The channels on the terrestrial delta top (i.e., p_land) move slightly slower in the treatment compared to the control. Lateral channel movement increases (i.e., channels move faster) slightly as the channels approach the mean shoreline in the treatment, whereas channel movement is less variable across the terrestrial delta top in the control. The perceived slower channel movement past the mean shoreline in both experiments is due to the fact that channels do not always reach that far into the basin (Fig. 4c), which artificially creates longer channel mobility timescales. This is why lateral mobility is typically analyzed on a basin-wide scale in experimental river deltas, as shown in Fig. 7a. Similar to the T_mob, the planform overlap metric (another proxy for an avulsion timescale; Wickert et al., 2013) shows little difference between control and treatment experiments (Appendix B; Fig. B5). However, it takes about one-half of the time to modify 90 % delta top by at least 1 mm of aggradation in the treatment compared to the control experiment (Fig. 7b), which agrees with the channel and far-field channel aggradation rates and is due to the addition of the non-fluvial proxy in the low-lying region of the delta top that induces far-field sedimentation in the treatment experiment that would otherwise not receive sediment, which is the case in the control experiment. Thus, the addition of non-fluvial sediment in the marsh window (∼ 0.25 mm h⁻¹ for areas in the stable marsh window) is more widespread and outpaces (on a short-term scale) the far-field sediment deposition of the control experiment. Relatedly, because the channels are the same depth (basin-wide mean; Eq. 1) and the experiments both have long-term, basin-wide aggradation rates of 0.25 mm h⁻¹ (or RSLR_b), the basin-wide T_c (Eq. 6) is the same (Table 2).

Figure 7 (a) The fraction of the delta not visited by a channel for the control (blue) and treatment (green) experiments versus time in run hours. (b) The fraction of the delta top that has not accumulated at least 1 mm of sediment for the control (blue) and treatment (green) experiments versus time in run hours. (c) The 50 % lateral channel mobility (h), or how long it takes the channels to visit 50 % of the radial transect, as a function of radial distance from the apex (m) for both experiments is shown. (d) The 90 % lateral channel mobility (h), or how long it takes the channels to visit 90 % of the radial transect, as a function of radial distance from the apex (m) for both experiments is shown.

4 Discussion

4.1 The impact of non-fluvial deposition on deltaic channel evolution

The impact of the non-fluvial deposition on deltas is the making of longer channels that must traverse a slope break and a large platform near sea level (Figs. 4 and 8). The channels in the control experiment are wider on average but typically turn into sheet flow before they reach the ocean basin (Fig. 8a; Appendix B; Fig. B3). Near the shoreline, the delta top is usually dominated by one steeply sloping channel that contains a large amount of overbank flow. As such, we observe topographic flow expansions (Sittoni et al., 2014) and distally increasing deposition rates (Fig. 6a), indicating a morphodynamic backwater control on channel filling and migration. Because of the high deposition rates in the lobes, the channels lose confinement before they reach the local shoreline, leading to abrupt backstepping. The flow expansions occur in the distal channel lobes but impact (i.e., decrease) the channelization of the entire system, leading to 60 % of the terrestrial delta being covered with flow on average. This mechanism of channel movement and fluvial sedimentation produces channels that only reach the local shoreline about 24 % of the time. The resulting overbank flow has channel beds that remain above sea level to the shoreline, producing a backwater length that is often effectively zero. Thus, topographic flow expansions are the primary control on channel kinematics in the control experiment, which is similar to previous physical experiments without variable discharge (e.g., Li et al., 2017; Hoyal and Sheets, 2009).

Figure 8 Control and treatment channel comparison. (a) Aerial image of the control delta top at run hour 81, with contour lines at 5 mm (green) relative to sea level (r.s.l.), sea level (black), and −9 mm r.s.l. (blue). (b) Aerial image of the treatment experiment at run hour 221, with contour lines at 5 mm r.s.l. (green), sea level (black), and −9 mm r.s.l. (blue). Note that the experiments are run in equilibrium, so both the control and treatment experiments shown here are in the same phase of evolution, despite being at different experimental run hours. We choose these two images to highlight the difference in the channel morphology.

A dynamic network of feedbacks created from the deposition of non-fluvial sediment allows the treatment experiment to exist in dynamic equilibrium with constant relative sea level rise but with a very different morphology than the control experiment. Near the shoreline, the treatment channels are gently sloping and contain less overbank flow than the control experiment. This leads to a significant decrease in channel deposition rate near the shoreline compared to the control (Fig. 6a). We observe a long hydrodynamic backwater reach during 95 % of the treatment experiment, despite constant water discharge (Fig. 5). Enhanced lobe progradation has been previously shown to increase backwater length (Moodie et al., 2019; Brooke et al., 2022), which agrees with our finding that non-fluvial deposition creates channels that extend further past the mean shoreline (creating elongated lobes) and increases the hydrodynamic backwater length. As such, the non-fluvial deposition, the reduced distal channel deposition rates, and the presence of a significant hydrodynamic backwater reach appear to be intimately related.

The treatment channels are also narrower (Fig. 4b), deeper (Fig. 6b), and longer (Fig. 4c) than the control channels. The treatment channels are narrower because the flow is partitioned between a larger number of distributaries. Though some flow is lost to overbank flooding in the distal portion of the delta (Fig. 3b), to conserve flow that becomes distributed over an increasing number of channels, the channels must become shallower and narrower with increasing distance from the apex (Figs. 4b, c and 6b). These findings agree with previous studies that show that deltas with vegetation and/or increased fine sediment fraction have increased channelization and narrower, longer, and deeper channels (Lauzon and Murray, 2018; Caldwell and Edmonds, 2014). While these studies primarily focus the hydrodynamic effects of vegetation-induced drag, we show that the extra mass accumulating in these vast wetland platforms impacts the channel morphology in similar ways. The formation of more distributary channels offsets the narrower channels in the treatment experiment to keep the channelized area similar between the two experiments. Furthermore, flow is more constrained to the channels (Fig. 3b), indicating channel aggradation as being the primary form of fluvial sedimentation.

Interestingly, the addition of non-fluvial sediment in the low-lying overbank regions of the treatment delta top creates channels that have a long backwater reach. It is often suggested that variable discharge is necessary to create such conditions (e.g., Lamb et al., 2012; Ganti et al., 2016a). Though we cannot rule out an influence of variable discharge in other settings, it was not the control on hydrodynamic backwater here. While these constant discharge experiments both show increases in channel aggradation rate and decreases in avulsion timescale with distance downstream, the treatment experiment lacks the significant increase in aggradation rate near the local shoreline that produced classic topographic flow expansions in the control experiment (Fig. 6a). These results suggest that the channel avulsions triggered by a peak in channel sedimentation occur less frequently in the treatment experiment compared to the control. Recent research shows that avulsion locations may occur because of geometric constraints (e.g., a break in slope) and can be determined without backwater hydrodynamics (Ratliff et al., 2021; Prasojo et al., 2022). Unlike the control, the treatment experiment exhibits a channel slope break near the local shoreline (Fig. 4a), and the channels usually extend past the local shoreline, creating significant backwater lengths often >  0.5 m (Fig. 8). Hence, we propose the following new control on the hydrodynamic backwater: non-fluvial sedimentation (i.e., wetland accretion). Thus, the hydrodynamic backwater promoted by non-fluvial deposition has a strong control on channel kinematics in the treatment experiment, which is similar to global deltas (e.g., Ratliff et al., 2021; Chatanantavet et al., 2012; Lamb et al., 2012). This has important implications for controls on avulsions and channel movement in field-scale deltas. Both experiments have temporal variability in backwater length that is roughly equivalent to the compensation and lateral channel movement timescales. Though outside the scope of this work, this is an area for further discovery.

Despite the widely recognized relationship between avulsion locations and backwater length, the experiments have vastly different backwater lengths and channel in-filling (avulsion) timescales (Table 2), with no clear difference in avulsion location or long-term channel mobility. From this, we can only conclude that mass balance and compensation are what control the long-term channel mobility here, as the mass balance (Sanks et al., 2022) and compensation (Table 2) were nearly identical between the experiments. The basin-wide timescale of lateral channel movement (T_mob) is similar in both experiments (Fig. 7a; Table 2). Previous research has shown that vegetation decreases the lateral mobility of channels, effectively slowing channel migration in experimental and field systems (Wickert et al., 2013; Tal and Paola, 2010; Ielpi and Lapôtre, 2020), due to increased channelization and enhanced bank stability. For this reason, we hypothesized that the treatment experiment would have a longer T_mob, but this is not the case. Two processes are likely at play here. First, the river sediment mixture includes a polymer for added cohesion, which was created to simulate the cohesive properties of vegetation (Hoyal and Sheets, 2009). As such, the proxy wetland (non-riverine) sediment does not add a considerable amount of cohesion to the system. Second, the volume of sediment that accumulates on the terrestrial delta top is similar in the two experiments (Sanks et al., 2022). Combined with the systems evolving under the same base level rise rates (Table 1), the long-term distribution of fluvial sediment needs to be similar. In other words, lateral mobility and channel movement is driven by the compensation timescale, so the systems cannot have long-term, basin-wide lateral mobility timescales that are much longer than the compensation timescale. Because the control experiment was stable in its channel dynamics (i.e., the lateral mobility timescale was longer than the compensation timescale), this limits the ability for the wetland sediment proxy to further decrease channel mobility. This is supported by the long-lasting (∼ 250 h) channel at the beginning of the treatment experiment that led to increased mobility in the remainder of the treatment experiment, which helped fill the accommodation and counteract the effects of the initial long-lasting channel (Fig. B8). Though we observe a smaller ratio of far-field aggradation to channel aggradation in the treatment experiment than the control experiment (by ∼ 50 %), and the maximum channel aggradation is significantly reduced, compaction and a set base level rise rate create two systems with a similar terrestrial mass balance over their lifespans (Sanks et al., 2022). Thus, it is actually unsurprising that the basin-wide T_mob is similar. If the proxy wetland material were added to an experimental delta, where the riverine material did not have added cohesion via a polymer, and if the lateral mobility timescale were considerably shorter than the compensation timescale, then the slower channel in-filling timescale in the treatment suggests that the proxy material may decrease the lateral mobility compared to a non-cohesive delta without proxy material. This is an avenue for further exploration. Despite these results, we observe subtle changes in lateral mobility as a function of distance from the apex (Fig. 7c and d), which shows that non-fluvial deposition in the treatment experiment decreases channel movement slightly for area of the delta that is less than 0.8 m from the apex compared to the control, leading to longer timescales necessary to reach 50 % and 90 % of this area. While interesting, more work is needed to fully understand lateral mobility as a function of distance in experimental-, numerical-, and field-scale river deltas.

We note again that our non-riverine sediment proxy is a simple proxy simulating the effect of mass that accumulates in river deltas through non-riverine processes in the coastal zone. We adapt a model tying primary production in wetlands to elevation relative to mean high tide (Morris et al., 2002) and deposit the proxy material in situ as a function of delta-top elevation relative to sea level. Though we neglect the effects of vegetation on hydrodynamics, vegetation has been shown to increase sediment trapping, slow down overbank flow, and create channel confinement (Nardin and Edmonds, 2014; Beltrán-Burgos et al., 2023). We observe similar results here without an explicit stem density. For example, the amount of overbank flow decreased in the treatment compared to the control, the channels extended further into the basin in the treatment compared to the control, and the channel in-filling rate decreased in the treatment compared to the control. While it is clear that ecological changes can profoundly affect the coastal zone, these results stress that the presence of marshes or tidal flats in the coastal zones of deltas exert a fundamental control on the delta's channel dynamics.

4.2 Implications

Planned river diversions are a critical management strategy for river deltas threatened by sea level rise. New river channels will be formed, and their morphology and kinematics will determine sediment deposition and accumulation; thus, understanding the interaction of fluvial and non-fluvial sedimentation is important for successful river diversions. While the distribution of fluvial sediment is extremely important on short timescales in these scenarios, the movement of the channels will control the long-term success of these engineering strategies. Our study suggests that channel and floodplain aggradation rates depend on the non-fluvial sedimentation occurring in the system. For example, wetland aggradation appears to decrease the rate of maximum channel aggradation and amount of overbank flow, so it is important to incorporate this into management plans.

We also show that non-fluvial sedimentation in the low-lying region of the delta top changes the morphology of the channels and is the primary control of hydrodynamic backwater effects in this experimental setting. Wetlands are deteriorating in many coastal deltas worldwide (e.g., Couvillion et al., 2017; Morton et al., 2015; Kirwan et al., 2010; Reed, 1995). If marshes are no longer present in the deltaic system, then it is possible that the morphology of the entire delta and its channels will change. Thus, we suggest that channel sediment deposition in a field-scale delta without marshes may be influenced more by topographic flow expansions and morphodynamic backwater effects than the hydrodynamic backwater effects that currently dominate most global deltas.

Furthermore, river avulsions are potentially catastrophic for communities (e.g., Slingerland and Smith, 2004), which may lead to the eventual regulation of hydrodynamics (e.g., flood stages) to control the avulsion location and river slope break location in deltaic systems, similar to engineered avulsions based on river sedimentation proposed by Moodie and Nittrouer (2021). We show that the marsh platform and its floodplain deposition produce the slope break and subsequent “hydrodynamic” backwater effect without unsteady hydrodynamics; thus, controlling the hydrodynamics of rivers may not be necessary.

5 Conclusions

The addition of the marsh proxy in the experimental setting produced some surprising dynamics. The channel morphology of the two systems varies greatly. We show that the in situ accumulation of sediment in wetland platforms and tidal flats increases the trunk channel depth and number of distributary channels in deltaic systems, which allows the channels to extend further into the basin. This non-fluvial deposition decreases the tendency for unchannelized (overbank) flow in the floodplain, effectively concentrating flow to the channels. In turn, the channels deposit sediment more consistently along the channel length, significantly decreasing the channel aggradation peak near the shoreline that occurs in a control experiment and is characteristic of morphodynamic backwater effects. Importantly, the presence of non-fluvial sedimentation in the treatment experiment creates a long hydrodynamic backwater reach without inducing variable discharge. Even though we observe these very clear differences in the channel morphology, the basin-wide timescale of lateral channel movement remains similar because it is fundamentally controlled by compensation and mass balance. The non-fluvial (marsh) sedimentation proxy produces a channel morphology of the treatment experiment more akin to that observed in global deltas, and thus the treatment experiment can be used to better understand the kinematics and long-term suitability of planned river diversions, avulsions, and general management strategies.

Appendix A: Channel maps

Channel morphology and kinematics provide key insight into the life cycle and fate of a river delta (e.g., Edmonds et al., 2009). To determine the channel morphology and kinematics in an experimental river delta, binary channel maps are needed. We use a software called ImageJ to hand map the channels in the experiments, as it retains basin coordinates. First, we load the RGB (red, green, and blue) image (collected every hour from the wet lidar scans) into the software (Fig. A1a) and then use the freehand selection tool to trace the outline of the channel map (Fig. A1b). Once the entire channel has been mapped, we export the channel map as a binary TIFF file (Edit → Selection → Create Mask). We use the resulting channel maps (Fig. A1c) to analyze the various properties and kinematics of both the control and treatment experiments. The control experiment channels were mapped by Ripul Dutt, and the treatment experiment channels were mapped by Kelly Sanks.

Figure A1 (a) An RGB image from hour 218 of the treatment experiment loaded into ImageJ. Note that the image size is 750×747 in the top-right corner, which shows the pixel dimensions of the lidar scans from the treatment experiment. (b) The yellow outline around the beginning of the channel illustrates one portion of channel map created using the freehand selection tool. The polygons were created by zooming into the channel and creating a polygon around the channel. (c) The resulting channel map overlain on the RGB image for hour 218 of the treatment experiment.

Appendix B: Channel kinematics

We use p_land as the delta area for most channel metric calculations, which describes the delta area that is above sea level for at least 50 % of the experiment. This area is smaller in the treatment experiment (∼ 2.0 m²) than the control (∼ 2.2 m²), but the size and locations at which the delta is above sea level in the treatment experiment are more variable (Fig. B1).

Figure B1 The area of land that is above sea level for various amounts of time for the experiments. The intersection of the black and blue line indicates the area of land (2.2 m²) that was above sea level for at least 50 % of the control experiment (p_land). Similarly, the intersection of the black and green line indicates p_land (2.0 m²) for the treatment experiment.

The lateral mobility and modification timescales come from the exponential decay function that describe how quickly the delta top (p_land) is visited by channels and how quickly the delta top accumulates sediment (Fig. B2). See Table 1 and Fig. 7a and b for the lateral mobility and modification in linear space.

Figure B2 (a) Fraction of the delta that has not been visited by a channel versus the measurement window in semilogy space (see Fig. 7a for linear space). (b) Fraction of the delta that has not accumulated at least 1 mm of sediment versus measurement window in semilogy space (see Fig. 7b for linear space).

We discuss two different channel movement styles in the results. These are the topographic flow expansions associated with a morphodynamic backwater and lateral movement associated with channel mouth extension and a hydrodynamic backwater. Topographic flow expansions are characterized by the mouth of the channel becoming unchannelized sheet flow over the entire depositional lobe (Sittoni et al., 2014). The control experiment is mostly characterized by topographic flow expansions, which lead to channel avulsions, likely because the channel beds are super elevated above relative sea level (Fig. B3). We see that the topographic flow expansions, which happen often in the control experiment and sometimes in the treatment experiment, are characterized by the formation of a main channel with some distributary flow. The distributary flow then begins to expand over the entire mouth bar (sheet flow), causing the main channel to recede, and the main channel begins to look for a new path at the final point of recession (avulsion location).

Figure B3 Topographic flow expansion from the control experiment displayed in hours 110, 113, 115, 117, 119, and 121. In this instance, it takes 11 h for the topographic flow expansion to occur.

Figure B4 Channel mouth bar formation and lateral switching from the treatment experiment displayed in hours 486, 489, 491, 496, 498, 500, 507, 515, and 520. In this example, it takes about 34 h for lateral movement of main channel.

While the treatment experiment does have topographic flow expansions, the channel movement is mostly characterized by the formation of long-lasting channel mouth bars and lateral switching of the channel paths (Fig. B4). We see the formation of a main channel, which occurred due to channel path finding and then lateral movement of the main channel to the new channel location. Upon the initiation of the main channel, it begins to extend further into the basin and past the mean shoreline and, in this example, stays in place for about 10 h. During this time, smaller channels form on the delta top in the “path-finding” phase. A small topographic flow expansion leads to the formation of two distributary channels over the mouth bar (not sheet flow like in the control), but ultimately, the real avulsion takes place near the entrance channel, where the main channel switches laterally, and the old channel is abandoned almost immediately. This process takes longer than the topographic flow expansion, but because of the path-finding phase, much of the delta is visited during this time (even though the main channel stays the same).

We look at the following channel kinematic metrics: lateral mobility (Li et al., 2017) and fraction of the delta unmodified (Li et al., 2017). We present the findings for the lateral mobility and fraction unmodified in Sect. 3. However, we describe three more channel movement metrics here. This includes the planform overlap metric (Wickert et al., 2013), the number of times a specific location (pixel) on the delta top was channelized, and the longest consecutive channelization for each pixel on the delta top.

We use the planform overlap as a check for the lateral mobility of the two systems (treatment and control). The planform overlap is similar to lateral mobility, as it is another commonly used proxy for the avulsion timescale in experimental settings (Wickert et al., 2013). The planform overlap is calculated by determining the time it takes for the channel maps to decorrelate from each other. To calculate this, we first calculate the number of changed pixels (channelized to unchannelized or unchannelized to channelized) for each time step (all 1 h time steps, all 2 h time steps, etc.). We see that there are fewer channel pixels changed in the treatment experiment (Fig. B5a). Following methods from Wickert et al. (2013), we scale the number of changed pixels by the number of pixels that would change via random scatter, which gives us ϕ (phi; Fig. B5b). We again see that there are fewer changed pixels in the treatment versus the control. We use this change to then determine essentially how long it takes for channel maps to decorrelate from each other. However, we see that the planform overlap tappers off at roughly the same time in both the control and treatment experiments. The planform overlap exhibits the same exponential decay trend shown by the lateral mobility metrics we used. Thus, we conclude that the lateral mobility of the two systems is similar (Fig. B5c).

Figure B5 Channel planform change. (a) The mean number of changed channel pixels (y) for various time steps (x) for the control (blue) and treatment (green). The standard deviation is shown as the light blue (control) and light green (treatment) polygons around the mean. (b) The mean number of changed pixels scaled by the number of changed pixels produced via random scatter for the control (blue) and treatment (green). The standard deviation about the mean is shown in light blue for the control and light green for the treatment. (c) The mean time it takes for the channel maps to decorrelate from each other for various time steps for both the control (blue) and treatment (green). The standard deviation about the mean is shown as the light blue (control) and light green (treatment) polygons.

While the lateral mobility timescales do not differ, we do see some subtle differences in channel mobility between the two experiments (Fig. B6). We see that there are more areas on the delta top in the treatment experiment that are rarely visited by a channel, which can be seen by the presence of more dark blue (less visited) areas on the delta top in Fig. B6b (treatment) than in Fig. B6a (control). This is supported by the probability distribution of the amount of time that each pixel was visited by a channel for the control and treatment experiments (Fig. B6c). The median time a pixel in the control experiment is visited by a channel is 37 h, but it is slightly shorter in the treatment experiment at 32 h. Furthermore, we see that more pixels are rarely visited by a channel in the treatment experiment than the control. This suggests that the treatment channels may have some amount of preferential flow (i.e., channels are more likely to reoccupy areas they have already visited), or that there is decreased lateral channel mobility combined with the efficient sweeping of the channels (supported by the longer channel in-filling timescale in the treatment), even though there is not a decrease in the time it takes the channels to move across the delta top over the delta's life cycle. We also observe a decrease in the channel mobility timescale in the treatment experiment due to the presence of a long-lasting channel at the beginning of the experiment (Fig. B8). It is possible that the long-lasting channel at the beginning of the treatment experiment impacted the average lateral mobility timescale, as the channels had to speed up to fill in the accommodation created by RSLR_b.

There is also no clear difference in the amount of consecutive time that channels occupied a specific area (Fig. B7a and b). The main distinction is that channels were able to consecutively occupy an area further from the entrance channel in the treatment experiment than the control experiment. When normalized by the total amount of time a channel occupied that pixel, we again show no significant difference. This normalization shows us that the channel tips tend to be only occupied once, whereas areas near the entrance tend to be occupied many different times (Fig. B7c and d). Finally, we see that a similar pattern emerges in the probability distribution of consecutive channelization (hours; Fig. B7e) compared to total channelization (hours) at a pixel (Fig. B6c). There is more area that is not visited by a channel often in the treatment experiment than the control but roughly the same amount of area visited by a channel for 50 consecutive hours or more (control = ∼ 3800 pixels and treatment = ∼ 3000 pixels; area not shown in Fig. B7e).

Figure B6 The time that each pixel in p_land is occupied by a channel. (a, b) The time that each pixel in the control (a) and treatment (b) experiments is occupied by a channel. (c) Probability distribution showing the time each pixel is visited by a channel for the control (blue) and treatment (green) experiments.

Figure B7 (a, b) The maximum amount of time each pixel in the (a) control and (b) treatment experiments was consecutively channelized (hours) for the entire channelized portion of the delta top. (c, d) The maximum amount of time each pixel in the control experiment was consecutively channelized (hours) divided by the total time each pixel was channelized in the (c) control and (d) treatment for the entire channelized portion of the delta top. (e) The probability distribution for the maximum consecutive time channelized (hours) for the control (blue) and treatment (green).

Last, we determine the time it takes to approach one e-folding reduction in non-channelized area from each run hour. The average fraction of the unchannelized area on the terrestrial delta in the control experiment is ∼ 0.91 ($\mathrm{1}-\frac{A_{\mathrm{c}}}{p_{\mathrm{land}}}$). To obtain the e-folding reduction, we divide the unchannelized fraction by e (2.71828), which is ∼ 0.34. Thus, we calculate the time it takes for channels to visit 66 % (1–0.34) of the control p_land. Similarly, the fraction of unchannelized area in the treatment experiment is ∼ 0.89. When reducing by an e-fold, we obtain 0.33, so we calculate the time it takes for channels to visit 67 % of the treatment p_land. We do this in order to normalize the channel visitation for differences in the channel and delta area in the control and treatment experiments. We see that the control experiment oscillates between about 50 and 100 h. Interestingly, the first ∼ 200 h of the treatment experiment have very slow channel mobility (e.g., taking between 80 and 150 h to visit 67 % of the p_land), but the last ∼ 300 h have very fast lateral channel mobility (e.g., taking between 20 and 60 h to visit 67 % of the p_land; Fig. B8). We see that the long-lasting channel(s) at the beginning of the treatment experiment significantly impacted the channel mobility in the second half of the experiment, potentially speeding up the channel motion, so that the delta top could stay in equilibrium.

Figure B8 The time to reach an e-folding reduction in a non-channelized area of the terrestrial delta in the control (blue) and treatment (green) from each run hour.

Appendix C: Figure 6

Below, we show Fig. 6 with all data. We note that control channel aggradation rate that is > 7 mm h⁻¹ (Fig. C1a) is likely a result of not many channels reaching that distance from the apex (Fig. 4c). This is also the case for the control channel in-filling rate that is > 7 mm h⁻¹ (Fig. C1c).

Figure C1 (a) The mean channel (solid line) and far-field (dashed line) aggradation rates (mm h⁻¹) in the control (blue) and treatment (green) experiments as a function of radial distance from the apex (m). (b) The mean trunk channel depth (mm) as a function of radial distance from the apex (m). (c) The aggradation difference (mm h⁻¹; i.e., the channel in-filling rate) between far-field and channel aggradation rates as a function of radial distance downstream (m). (d) The channel in-filling timescale (h) as a function of radial distance downstream (m).