Ancient Channel-Mouth Bifurcation Angles on Earth and Mars

INTRODUCTION

Understanding of the morphodynamic controls and topological complexities of distributary channel networks on river deltas has advanced rapidly over the past decade (e.g., Edmonds, Paola, et al., 2011; Hariharan et al., 2022; Konkol et al., 2022; Piliouras & Rowland, 2020; Tejedor et al., 2017). However, the translation of these insights from modern networks that are readily imaged via remote sensing to ancient ones remains difficult. The principal difficulty in this translation is the requirement of a map of network shape over length-scales of many channel-widths, with confidence that the shape corresponds to a geologic instant rather than an amalgamation of time. While this is rare with outcrop-based studies (but see Li & Bhattacharya, 2014), network shapes can be resolved in both 3D seismic volumes (e.g., Martin et al., 2018; Zeng et al., 2013) and in locations where preferential erosion of fine-grained floodplain sediments has exposed coarse-grained channel networks (Pain et al., 2007; Pain & Oilier, 1995), such as at various locations on Mars observed through remotely sensed/orbital imagery (e.g., DiBiase et al., 2013; Hughes et al., 2019; Malin & Edgett, 2003; Moore et al., 2003). These datasets allow for quantitative measurement that adds an important deep-time and planetary perspective on channel network dynamics. Here, we conduct such an analysis, combining terrestrial 3D seismic with Mars orbital imagery to measure and characterize a distribution of ancient channel mouth bifurcations, where a single channel tip split into two distributary channels during progradation (Edmonds & Slingerland, 2007; Fidolini & Ghinassi, 2016; Shaw & Mohrig, 2014; Wright, 1977).

Channel mouth bifurcations are a special case where a local network phenomenon offers promising clues to formative dynamics. Under the assumptions of a thin channel supplying flow to a semi-infinite basin, a channel bifurcation will grow toward an equilibrium angle of 72° to maximize flux at the two channel tips (Devauchelle et al., 2012; Ke et al., 2019). Despite standard errors larger than 15° for all datasets, empirical measurements have validated this observation in seepage tributary channels in Florida (71.9° ± 0.9° mean plus 95% confidence interval; n=4,966, Devauchelle et al., 2012), tributary networks forming in humid, low slope conditions across Earth (Seybold et al., 2017, 2018), in distributary networks found on field scale river deltas (70.4° ± 2.6°; n=222, Coffey & Shaw, 2017), and in experimental river deltas (68.3° ± 8.7°; n=21, Coffey & Shaw, 2017 and 74.1° ± 7.7°; n=13, Federici & Paola, 2003).

While the tendency toward this critical angle is remarkably consistent among channel mouth bifurcations, it is by no means universal. Angles tend to be smaller for bifurcations around braid bars and at channel avulsions. Dong et al. (2020) found a mean bifurcation angles of 47° (17° standard error) on 12 bars in the braided reach of the Selenga River Delta, Russia. Burge (2006) measured bifurcations on seven avulsion channels in gravel bedded rivers in New Brunswick, CA and found a mean angle of 48° (8° standard error). In numerical models of channel bifurcation, the width of channels, critical shear stress of sediment motion, and sediment transport nonlinearity used can influence the emergent bifurcation angle and other network properties (Ke et al., 2019; Żukowski et al., 2022). Even assuming the thin-finger theory responsible for the critical bifurcation angle of 72°, models that optimize growth rate of channels rather than even flow distribution at channel tips suggest a bifurcation angle of 78.5° (Devauchelle et al., 2017). Finally, the stability of bifurcation angles during network abandonment has not been investigated systematically, but significant angle widening directly prior to abandonment has been documented in particular cases (Cole et al., 2021; van Heerden & Roberts, 1988). Hence, there is no guarantee that ancient channel bifurcations will cluster at 72°.

The stratigraphic record of Earth is awash with distributary channel deposits (e.g. Fisher & McGowen, 1967; Korus & Joeckel, 2023; Payenberg et al., 2003, 2024; Zeng et al., 2001, 2013; see Figure 1a). Ancient deltaic deposits have also been interpreted in many locations on Mars (e.g., Burr et al., 2009; De Toffoli et al., 2021) and specifically for this study, in Aeolis Dorsa (DiBiase et al., 2013; Hughes et al., 2019) and Eberswalde crater (Malin & Edgett, 2003; Moore et al., 2003). Preservation of individual channel mouth bifurcations, however, is less commonplace and interpretation of those angles in comparison to other divergent features such as avulsions is challenging. If predictable bifurcation angles are preserved in the stratigraphic record, they could provide a potential tool for improving stratigraphic interpretation and forward modeling.

Figure 1.Stratigraphic distributary channel networks in (a) seismic coherence horizon slice from Pleistocene paleo-Mississippi delta deposits near South Marsh Island, Louisiana; and (b) stereo-derived color digital elevation model overlain on HiRISE imagery of martian topographically-inverted channel deposits preserved in Eberswalde crater (McEwen et al., 2007). Examples of interpreted amalgamated belts and multi-story belts (which were excluded from analysis) and individual channels (the focus of this analysis) are indicated.

METHODS

We measured the angle of interpreted channel mouth bifurcations preserved in terrestrial seismic volumes and imagery of Mars stratigraphy. The identification and interpretation of channel mouth bifurcations is a non-trivial task; the distinction between channel-mouth bifurcations, intersecting channel bodies, and channel avulsions is not easily made in the stratigraphic record (Li & Bhattacharya, 2014; Olariu & Bhattacharya, 2006). We relied on three criteria for mapping. First, we required that the channel deposits were not observed to rejoin downstream to prevent the measurement of braiding channels (Figure 1B). Second, channel deposits must be found in the same stratigraphic story to limit the likelihood of avulsion or reoccupation. Third, we required that the measured channel bodies were necessarily the result of aggradation within a single channel, rather than the amalgamation of many channels into a channel belt. This criterion allows us to directly estimate channel width from deposit width.

The inverted channel features of the Cedar Mountain (Utah, USA) and Guadalope-Mantarranya (Spain) Formations are used as analogs for many topographically inverted sedimentary features on Mars, and have been shown to be the amalgamation of many channels (e.g., Cardenas et al., 2020; Hayden et al., 2021). The variable belt width and internal scroll bar morphology of proximal delta deposits (lower left of figure 1B) also suggest amalgamation and are not analyzed here. However, distal deposits on Mars (figure 2B) show narrow channel bodies with small width variation, suggesting a single channel that was not reoccupied after abandonment. Indeed, many terrestrial rivers with amalgamated sandy channel belts narrow near the coastline such that the sand deposit width and channel width are subequal (Fernandes et al., 2016; Martin et al., 2018; Wahbi & Blum, 2023). Ongoing investigation of distal channel bodies reveals that their internal stratigraphy can be consistent with single channel occupation and negligible lateral migration (Hughes et al., 2024). These criteria give us high confidence that the features we are interpreting are indeed ancient channel mouth bifurcations.

Figure 2.Example of determination of uncertainty in angle measurements. Panel a shows sweetness measured at Marsh Island Louisiana. Panel b shows hillshade overlain digital elevation model at Aeolis Dorsa, Mars. Black lines show maximum angle and red lines show minimum angle measurements.

Data for distributary networks found in seismic volumes from Earth were compiled from published sources of 3D seismic data from distributary channel strata (Y. Dong et al., 2017; Hao et al., 2014; Hart, 2008; Zeng et al., 2001, 2012, 2013), as well as one set of 3D seismic data from South Marsh Island (B-31-92-LA), northern Gulf of Mexico that was analyzed by the authors. This survey covers an area of 1,196 km² and was collected in 1992 by Bureau of Ocean Energy Management. We searched for channel mouth bifurcations across a significant fraction of all publicly available seismic volumes from the Gulf of Mexico and New Zealand (Table 1). Despite the prevalence of fluvial channel bodies within these systems (Alqahtani et al., 2015; Armstrong et al., 2014; Cardenas et al., 2022), bifurcations matching the criteria above were rare. The strata where channel-mouth bifurcations were found are interpreted predominantly as shallow prograding river delta systems in both marine and lacustrine basins. 109 bifurcations were measured from 9 separate seismic data sets. Due to the proprietary nature of many of the seismic data sets, angles had to be measured from published images, rather than the seismic datasets themselves for all cases except South Marsh Island.

Interpreted channel-mouth bifurcations were compiled for martian deltas using publicly available imagery and stereo-derived digital terrain models (DTMs) and orthophotos collected by the High Resolution Imaging Science Experiment (HiRISE) and the Context Camera (CTX) on the Mars Reconnaissance Orbiter (Observations IDs listed in Table 1; Malin et al., 2007; McEwen et al., 2007). Mars DTMs are products of the NASA Ames Stereo Pipeline (ASP) (Beyer et al., 2018; Shean et al., 2016) where stereo-pairs were tied to Mars Orbiter Laser Altimeter (MOLA) point shot data to minimize error associated with regional slope. HiRISE DTMs have pixel dimensions of ~1 meter per side, and orthophotos are ~0.25 meters per side. CTX DTMs have pixel dimensions of ~18 meters per side, and orthophotos are ~6 meters per side. Deltas with recognizable distributary networks previously identified at Eberswalde crater (see Wood, 2006) and in the Aeolis Dorsa region (DiBiase et al., 2013; Hughes et al., 2019) were used in our analysis. 50 bifurcations were measured from these three localities.

Angles were measured using the methods of Coffey and Shaw (2017): by selecting the bifurcation apex and margins of the two daughter channels along the inside of the bifurcation (see Figure 2). This method is important, because channel mouth bifurcation angles tend to decrease over longer measurement scales but are consistent with 72° for length scales less than 1 upstream channel width. Upstream channel width was not typically apparent in stratigraphic data, and visible channel margins in the downstream were often incomplete, limiting our ability to measure at locations 1 parent channel width downstream from the apex; however, we did our best to measure over an approximately similar scale to that of the channel widths.

To account for the difficulty with estimating channel boundaries and other sources of uncertainty, we found it appropriate to map an envelope of possible bifurcation angle. Hence, we measured a minimum and maximum possible angle at each site. In terrestrial seismic data, uncertainties were associated with unclear transitions in the time or horizon slice of an attribute between channels and islands. On martian outcrop data, channel boundaries were mapped using the slope break of singular ridges. Mapping the apex of the bifurcation was particularly difficult in this case. In Mars outcrop data, uncertainties were mostly associated with identifying the exact interior margins of bifurcate channels (Figure 2).

We compared angle distributions using the estimated mean and associated 95% confidence interval of channel mouth bifurcation angle preserved in stratigraphy to the theoretically predicted bifurcation angle (72°). Monte Carlo simulations were conducted to account for both a limited number of samples and uncertainty in each measurement (Rubinstein & Kroese, 2016). Sampled distributions were bootstrapped with replacement to produce 10⁶ sets of n angle envelopes (n = 109 and 50 for terrestrial and martian datasets, respectively). For each set, a simulated mean was calculated by randomly selecting angles from each envelope assuming a uniform distribution between the bounds. While the output of these simulations produces a mean angle which converges to the mean of the midpoints of all measurement envelopes, the confidence interval is related to the spread of the envelopes. Matlab code for the implementation of the Monte Carlo simulations is provided in the Supplemental Materials.

RESULTS

Comparison between modern terrestrial channel mouth bifurcations from Coffey and Shaw (2017), ancient terrestrial bifurcations derived from seismic volumes, and ancient Martian bifurcations derived from remotely sensed outcrops suggests similarity between terrestrial datasets, but slightly larger angles on Mars (Figures 3 and 4).

The 109 bifurcation angles measured from terrestrial seismic data had an associated envelope size range of 1.7° to 58.3° (mean 26.6°; Figure 3a). Monte Carlo estimation of the mean angle yielded 71.9° ± 3.8° (95% confidence interval on mean) with an expected standard error of 18.5° (Figure 4d). A rank-sum U test (Mann & Whitney, 1947; Wilcoxon, 1992) on the mean values for each bifurcation (the midpoint between the minimum and maximum angle interpretation) could not reject the hypothesis that the samples from terrestrial seismic and modern terrestrial deltas (Figure 4a, b) were derived from the same distribution (p-value = 0.40). This is because the standard error of modern Earth angles (18.6°) is similar to those found in terrestrial seismic.

Figure 3.(a) Each line represents the measurement of a bifurcation, connecting the minimum and maximum angles. Terrestrial seismic measurements are in gray and martian measurements are in red. (b) p(angle) shows the fraction of all measurements for which the angle falls between the minimum and maximum mapped angle.

Figure 4.Histograms of bifurcation angles derived from (a) modern and experimental deltas (after Coffey & Shaw, 2017), and from midpoints of minimum and maximum interpreted bifurcation angles from stratigraphic sequences derived from (b) terrestrial seismic data and (c) Martian orbital imagery. (d) and (e) show distributions of sample means from bootstrapped, Monte Carlo (MC) simulations from Earth and Mars, respectively.

The 50 bifurcation angles measured from martian outcrop data had an associated uncertainty ranging from 2.4° to 41.2° (mean 18.2°; Figure 3a). Monte Carlo estimation of the mean angle yielded 80.1° ± 4.8° (95% confidence interval on mean) with an expected standard error of 16.4° (Figure 4e). This prediction is 11% larger than for the terrestrial seismic data. A rank-sum U test on the mean values for each bifurcation rejected the hypothesis that the samples from martian outcrop and modern terrestrial deltas were derived from the same distribution (p-value = 0.011).

The fraction of measurements consistent with any given angle was compared between terrestrial seismic and Mars outcrop (Figure 3b). The distribution of maximum angles is similar between Mars outcrop and terrestrial seismic: in both datasets, a bifurcation angle of 84° is consistent with about 40% of measurements and an angle of 108° is consistent with about 10% of measurements. In contrast, larger minimum angles were mapped on Mars outcrops: a bifurcation angle of 60° was consistent with 43% of measurements in terrestrial seismic and 22° or martian outcrop measurements. A bifurcation angle of 36° was consistent with 9% of measurements in terrestrial seismic and was never measured on martial outcrops (0% consistent).

In summary, we find no significant difference between the channel mouth bifurcation angles on modern Earth and terrestrial seismic data. However, measurements from martian data deviate from prediction and from the distributions found in modern deltas and terrestrial seismic data.

DISCUSSION

Distributions of bifurcation angles measured from terrestrial seismic volumes are found to be statistically similar to modern river delta channel mouth bifurcations (Coffey & Shaw, 2017). Further, these mean values are consistent with theory for bifurcation initiation of thin channels in a diffusive flow field (Devauchelle et al., 2012). From this we can infer that distributary channels are remarkably stable through their abandonment. On Earth at least, relict channel mouth bifurcations preserve the dynamics that formed them with quantitative reproducibility despite their scarcity and the challenges inherent in interpretation and measurement.

Bifurcation angles measured from martian outcrops were statistically wider than Earth data. This difference can be explained either by measurement bias between outcrop and seismic data, or by differences in sediment transport between the planets. Sampling bias could exist given that the measured bifurcations are erosional remnants of channel deposits, and that the erosional process would tend toward widening of these angles as ridge edges erode back in a manner that increased the minimum mapped bifurcation angle (Figure 3b), yielding larger estimated mean angle values. Similar maximum and larger minimum angles on martian outcrop relative to terrestrial seismic data (Figure 3) supports this interpretation. In this case, an additional model reconstructing the interior channel margin is required to estimate true channel bifurcation.

On the other hand, the planetary environment of Mars could conceivably influence distributary channel bifurcation dynamics. The critical bifurcation angle of 72° derived by Devauchelle et al. (2012) does not depend on gravity, fluid density, or sediment transport rates, but channels must be thin (i.e. their width is so small relative to their length that they can be approximated by a single line). Models that investigate network growth with channels of finite width show that bifurcation angle is sensitive to critical Shields stress $\tau_{cr}^{*}$ and transport nonlinearity $\alpha$ that compose a general bed-material transport equation (Fernandez Luque & Van Beek, 1976; García, 2008).

$$q_{s} = CD\left( \tau_{b}^{*} - \tau_{cr}^{*} \right)^{\alpha}\tag{1}$$

Here, q_s is sediment flux per unit width, C is a coefficient related to resistance, D is grain diameter, $\tau_{b}^{*}$ is dimensionless shear stress ($\tau_{b}^{*} = \tau_{b}/(\rho_{s} - \rho)gD$ with g acceleration due to gravity and $\rho_{s}$ and $\rho$ are the sediment and fluid densities). Numerical models of prograding distributary channel networks by Ke et al. (2019) showed that values of $\alpha = 1.5$ and $\tau_{cr}^{*} = 0.1$ recreated network characteristics of the prograding, bifurcating, Wax Lake Delta, Louisiana with a numerically modeled mean bifurcation angle of 68.7°, similar to direct angle measurements for that particular system (mean 70.8°; Coffey & Shaw, 2017). However, when $\alpha$ and $\tau_{cr}^{*}$ were increased, the mean bifurcation angle also increased (Figure 5, Appendix). Holding $\alpha = 1.5$ (the standard value for bedload transport on Earth), an increase of $\tau_{cr}^{*}$ from 0.10 to 0.31 would increase the expected bifurcation angle by the 11% between Earth and Mars. Notably, $\tau_{cr}^{*}$ can be quite variable even on Earth (Otsubo & Muraoka, 1988; Parchure & Mehta, 1985), but a fundamental change in sediment erodibility on Mars could alter the bifurcation angle.

Figure 5.Predictions of mean bifurcation angle (color field and contours) derived from numerical models of Ke et al. (2019) as a function of $\alpha$ and $\tau_{cr}^{*}$. The gray and red zones show the 95% confidence interval of mean bifurcation angle on Earth and Mars.

Holding $\tau_{cr}^{*}$ constant at 0.1, the 11% larger mean angle on Mars could also be explained by increasing $\alpha$ from 1.50 to 1.73. On Earth, $\alpha$ is well constrained for sandy bed material on Earth at 1.5-1.6 (Wong & Parker, 2006), but reaches 1.68 for the Huanghe River (China) due to its very fine sand bed-material transported in suspension and lack of large in-channel bedforms (Ma et al., 2017). The reduced gravity on Mars is expected to increase the amount of bed-material transport in suspension, particularly for grainsizes near the threshold of suspension (Braat et al., 2024), conceivably yielding an $\alpha$ value that is consistently larger on Mars. While channel characteristics like width-to-depth ratio and water surface slope are expected to change between Earth and Mars, it is the sediment transport nonlinearity terms in Eq. (1) that functionally influence the emergent bifurcation angle and network growth laws (Ke et al., 2019; Żukowski et al., 2022). Data do not yet exist to test the influence of $\alpha$ and $\tau_{cr}^{*}$ on modern Earth deltas, but it is conceivable to test this influence taking advantage of these edge cases on Earth. Hence, while we consider measurement bias from outcrops to be a plausible explanation, increased bifurcation angles from altered sediment transport processes on other planets cannot be ruled out.

This analysis straddles the gulf between modern and ancient datasets. In the case of terrestrial seismic data, our results show that distributary network properties can be mapped to quantitative accuracy in the subsurface. While the sample size of the terrestrial seismic and martian outcrop datasets will likely increase with the availability of seismic volumes and martian digital elevations models, we do not anticipate these numbers will ever approach the number of angles that can be measured on a single remote sensing platform of modern Earth. Additionally, uncertainty in bifurcation angle measurement will always be larger when mapped in ancient strata. Finally, the lack of widespread channel mouth bifurcations preserved in inverted relief on Earth or the present lack of subsurface basin imaging on Mars prevent a planetary comparison with controlled methods, although numerical models could conceivably address this problem. Despite these limitations, the results presented here provide a path forward for interpreting channel mouth bifurcations in ancient strata.

We expected that channel mouth bifurcations would be a favorable environment for stratigraphic preservation; yet bifurcations were remarkably rare in seismic volumes. We speculate on possible reasons why they are so rare. First, the low accommodation, progradational, environments where channel mouth bifurcations are found are liable to be eroded immediately after deposition by a prograding fluvial system. Second, seismic reflection requires a contrast between the channel and inter-channel facies, so sandy bifurcating networks in a sandy system may not be resolvable. Third, lacustrine basins are overrepresented in our database (Table 1), suggesting that some combination of shallow basin depths, lateral confinement, or protection from basin processes may play an outsized role in bifurcation growth and preservation.

Implications for subsurface mapping and martian exploration

Understanding of the persistence of this morphology through channel network evolution, abandonment, preservation, and limited measurement sample sizes suggests it can be used to predict the stratigraphic architecture of sedimentary reservoirs. For example, a distribution of bifurcation angles may be a key variable in a training image (Scheidt et al., 2016), or in a process-mimicking statistical model (Pyrcz & Deutsch, 2014). Additionally, the persistence of this morphology is promising in its potential to be distinguished from other types of branching networks. While the review of avulsion angles is not systematic or large, a mean of 47° found in our literature review appears to be significantly smaller than the channel mouth bifurcations studied here.

CONCLUSION

By quantitatively mapping ancient channel mouth bifurcations on Earth and Mars, we were able to test models of channel network dynamics in ancient environments. Our results indicate that the angles of river delta channel mouth bifurcations interpreted in the terrestrial stratigraphic record are consistent with those observed on modern river deltas. This angle, with a mean of 71.9°, is consistent with that predicted from theory of network growth in a diffusive flow field. Bifurcations on martian outcrops have wider angles, with a mean of 80.1°. This difference could either be from measurement bias introduced by outcrops, or a slightly larger bifurcation angle due to differing sediment transport nonlinearity on Mars. Ultimately, these findings yield insight into the processes operating at river mouth bifurcations, as well as a predictive tool for stratigraphic modeling.

ACKNOWLEDGMENTS

We are grateful for comments from Edwin Kite and Tim Goudge who provided both constructive reviews and digital terrain models from Mars orbital imagery as well as for comments from an anonymous reviewer. We thank Heather Bedle for helpful discussions. This work was supported by DOE grant DESC0016163 to J.S. R.M. was partially supported through an NCED2 postdoctoral fellowship (through NSF-EAR-1246761), and funds from the Maurice F. Storm Endowed Chair of Petroleum Geology provided by Christopher Liner.