Chaos in a simple model of a delta network.

The flux partitioning in delta networks controls how deltas build land and generate stratigraphy. Here, we study flux-partitioning dynamics in a delta network using a simple numerical model consisting of two orders of bifurcations. Previous work on single bifurcations has shown periodic behavior arising due to the interplay between channel deepening and downstream deposition. We find that coupling between upstream and downstream bifurcations can lead to chaos; despite its simplicity, our model generates surprisingly complex aperiodic yet bounded dynamics. Our model exhibits sensitive dependence on initial conditions, the hallmark signature of chaos, implying long-term unpredictability of delta networks. However, estimates of the predictability horizon suggest substantial room for improvement in delta-network modeling before fundamental limits on predictability are encountered. We also observe periodic windows, implying that a change in forcing (e.g., due to climate change) could cause a delta to switch from predictable to unpredictable or vice versa. We test our model by using it to generate stratigraphy; converting the temporal Lyapunov exponent to vertical distance using the mean sedimentation rate, we observe qualitatively realistic patterns such as upwards fining and scale-dependent compensation statistics, consistent with ancient and experimental systems. We suggest that chaotic behavior may be common in geomorphic systems and that it implies fundamental bounds on their predictability. We conclude that while delta "weather" (precise configuration) is unpredictable in the long-term, delta "climate" (statistical behavior) is predictable.

Deltas are landforms arising from the deposition of sediment by a river entering a standing body of water. Deltas worldwide are highly populated and widely relied upon for agriculture and navigation but are under threat of collapse due to relative sea-level rise and sediment-supply reduction (1–3). Additionally, delta deposits are important reservoirs for groundwater and hydrocarbons, and better understanding their stratigraphic architecture would improve prediction of subsurface fluid flow (4).

Delta networks deliver sediment to different parts of the delta, controlling where land is built. Deltas have been proposed to self-organize their flux partitioning for a given network topology (5). The flux partitioning in the network can change over time through the process of avulsion, which may result in the abandonment of old channels and creation of new ones or merely a shift in the flux partitioning without modification of the number of channels (termed “soft avulsion”) (6). A wide body of research has aimed to better understand avulsions in order to better predict their timing and location and to explore how deltas record themselves in stratigraphy (7–14).

Our ability to accurately predict the future evolution of river deltas is hampered by our incomplete representation of the physical processes in our models, as well as uncertainty in boundary conditions. A similar problem arises in the prediction of subsurface stratigraphy based on incomplete data. Is there a fundamental limit to our ability to accurately predict delta evolution, which would apply even if we could perfectly simulate the relevant physical processes and boundary conditions?

Complex and irregular behavior, typically assumed to be stochastic, can be found everywhere in geomorphology. Examples include the motions of bedload particles (15), the evolution of dune fields (16), the organization of drainage networks (17, 18), the dynamics of braided (19) and meandering (20) rivers, and the kinematics of delta surfaces (21). Reproducing in detail the irregular dynamics of geomorphic systems is beyond the reach of existing models, but it is unknown to what extent detailed prediction of these systems is impossible rather than “merely” difficult.

In chaotic systems, long-term prediction is fundamentally impossible. Strogatz (22) proposed the following definition of chaos: “Chaos is aperiodic long-term behavior in a deterministic system that exhibits sensitive dependence on initial conditions.” Sensitive dependence on initial conditions refers to the idea that two systems started with nearly identical initial conditions diverge exponentially. This means that any error in our measurement of a system’s initial condition will grow exponentially through time. Because we can never measure a system with infinite precision, this implies a fundamental limit to long-term prediction. In contrast to stochastic systems, unpredictability arises in chaotic systems despite their lack of randomness; chaotic systems are by definition deterministic.

While many geomorphic systems display highly complex dynamics, demonstrating chaotic behavior is difficult. Major challenges include the lack of very long time series, which are typically required to detect chaos (23–25), and the fact that some methods to detect chaos can be fooled by certain types of noise (26). In the study of geomorphic systems, some have attempted to use spatial information to supplement the lack of temporal information needed to detect chaos (27), for example, from a time series of ripple migration (28).

Despite the complex behavior of natural geomorphic systems, a growing body of research on so-called “reduced-complexity models” has shown that complex dynamics that are at least qualitatively similar to those of natural systems can be generated even by models obeying relatively simple rules (16, 19, 29–34). Demonstration of chaotic dynamics in an idealized and simplified model of a natural system can be a strong suggestion that the system is chaotic; for example, the Lorenz equations are an extreme simplification of atmospheric convection but are believed to indicate that the atmosphere is itself chaotic (35). Chaotic behavior has been uncovered in models of geological phenomena such as earthquakes (36) and the geodynamo (37). Dynamic stratigraphic models are typically nonlinear dissipative systems and, therefore, are potentially susceptible to exhibiting chaos (38). More generally, the same could be said of many geomorphological models. However, to date, demonstrations of chaos in geomorphological models are rare: one example is from Pelletier (39), who identified deterministic chaos in a model of landform evolution, albeit with a definition of chaos that differs from the one we use in this paper. Another example is Phillips (40), who found chaos in a difference equation modeling regolith cover on a hillslope. Chaos has also been proposed for numerical models of river meandering (41, 42).

In this paper, we consider a numerical model of a simple delta network. The network consists of three coupled bifurcations: a channel splits into two branches at a bifurcation, each of which split into two additional branches. Our model is a further extension of the quasi-one dimensional (1D) bifurcation model developed by Bolla Pittaluga et al. (43) (hereafter, BRT) and extended by Salter et al. (44, 45) (hereafter, SPV) to include the effect of deposition in the downstream branches. We find that the coupling between the upstream and downstream bifurcations leads to chaotic dynamics in the flux partitioning of the network. This suggests that delta networks that have coupled bifurcations behave chaotically and that it is possible to quantitatively estimate the time scale over which long-term detailed prediction of their avulsion dynamics becomes impossible.

Model Overview

Our starting point is the BRT (43) quasi-1D bifurcation model, which has been widely used and built upon (46–52). This model allows for a transverse bed slope immediately upstream of the bifurcation, controlling the steering of water and sediment to the downstream branches. This is implemented by dividing the reach immediately upstream of the bifurcation into two laterally interacting cells, which have a length of , where is an order-one parameter and is the width of the upstream channel. The evolution of the bifurcation depends on the interplay between a positive and negative feedback: on one hand, if one of the downstream branches becomes deeper than the other, then in order to maintain a constant water level at the bifurcation, a larger discharge is required into the deeper branch. This exerts a nonlinear influence on the sediment transport rate, causing the deeper branch to become even deeper. Counteracting this positive feedback is the transverse sediment flux resulting from the presence of a cross-stream current and transverse bed slope immediately upstream of the bifurcation. Consistent with other existing bifurcation models (47, 53), the long-term behavior of the BRT model is always a fixed point, i.e., an unchanging state. A symmetric discharge partitioning is always a fixed point of the model, but depending on the interplay between the positive and negative feedback, it may be stable or unstable (51, 54). In the case that the symmetric fixed point is unstable, there exist two reciprocal asymmetric stable fixed points. This implies that following an arbitrarily small initial perturbation, an initially symmetric bifurcation tends to evolve toward a stable asymmetric configuration.

Whereas prior studies chose downstream boundary conditions resulting in sediment bypass (i.e., sediment output equal to sediment input), river deltas are landforms resulting from sediment deposition. SPV (44) imposed deposition via a specified bypass fraction , which is the fraction of the supplied sediment that exits the downstream end of the network. A bypass fraction implies deposition. Whereas bypass (F = 1) results in a fixed point, i.e., unchanging discharge partitioning, deposition creates the possibility of ongoing avulsion dynamics. The basic mechanism is as follows: because more sediment is routed to the dominant branch than the subordinate branch, its slope tends to become lower relative to the subordinate branch over time. Eventually, this slope difference causes a switch in the discharge partitioning. As with bypass, a symmetric fixed point always exists but may be stable or unstable. However, when it is unstable, rather than evolving toward a stable asymmetric fixed point, the bifurcation undergoes repeating avulsion dynamics, either soft avulsion (6) or full avulsion. These avulsion cycles are stable limit cycles, meaning that for a range of initial conditions, the bifurcation tends toward a periodic solution that is independent of the initial conditions.

To summarize, when the evolution of the bifurcation is decoupled from the downstream slopes, the system evolves toward a stable fixed point (43). In contrast, the coupling between the positive feedback in the bifurcation and deposition in the downstream branches results in stable limit cycles (44).

Here, we consider the scenario in which bifurcations are coupled in a simple network, with an upstream bifurcation coupled to bifurcations in each of its two downstream branches, as shown in Fig. 1. The model consists of a total of seven branches coupled via three bifurcations. The widths () and lengths () of the branches are fixed during model runs. While each bifurcation in isolation would produce a limit cycle, we find that the coupling between bifurcations can produce chaos. This is qualitatively analogous to the behavior of slider-block models of earthquakes: whereas a single slider-block is periodic, coupled slider-blocks introduce the possibility of chaotic dynamics (55).

Fig. 1.

(A) Buha River Delta, Qinghai Lake, China, N E. Note the presence of multiple orders of bifurcations and that the distributary branch widths are unequal, suggesting an asymmetric partitioning of fluxes at bifurcations. Image credit: Google Earth, © 2020 CNES/Airbus. (B) Schematic of the simple bifurcation-network model, which consists of a total of seven branches, coupled via three bifurcations as shown in the figure.

For simplicity, we assume that the geometry of the network (as described by and ) is symmetric about each of the three bifurcations. We note that relaxing this assumption does not qualitatively change the dynamics. We can cast our system in dimensionless terms by normalizing the network geometry ( and ) by the upstream width , normalizing bed elevation by the upstream channel depth , and normalizing fluxes by their value at the upstream boundary of the network. We introduce the system-wide mean deposition rate , computed as:where is the volumetric sediment discharge, is the bed porosity, and is the area of the delta-network topset, given in Eq. . This allows us to calculate the characteristic avulsion time scale, which we define as (56). We use this time scale to nondimensionalize time, i.e., . The dimensionless parameters required to specify run conditions are the normalized network geometry; the bypass fraction (Eq. ), which, in this paper, we will confine to in order to obtain net deposition; the upstream channel width-to-depth ratio ; Shields stress ; slope ; the exponent and critical Shield stress in a generic sediment transport formula (Eq. ); and the model parameters and , which specify the length of the divided reach and strength of bed load steering from a lateral bed slope, respectively. The parameter values used in this paper are listed in Table 1.

Table 1.

Model run parameters for Figs. 2–4

Figure	Fig. 2	Figs. 3 and 4
B_1/h_1	See figure	100
τ_*1	0.08	0.06
S_1	10^−3	10^−3
F	0	0.7
L_1/B_1	18	18
L_2/B_4	30	30
L_4/B_4	240	66
α	3	3
r	0.5	0.5
m	1.5	1.5
τ_*cr	0.047	0.047

Fig. 2.

Trajectories through a subset of state space, using water-discharge fraction through branches 4, 5, and 6 (color) as independent state variables. All model parameters are held constant except for the upstream width-to-depth ratio, which differs among A–E. (A) Example of a stable limit cycle (i.e., periodic solution). Notice that the cycle is not symmetric across the line . Inset shows the corresponding time series of discharge fraction through a single branch. (B) Example of a quasiperiodic solution. (C) Example of a stable limit cycle for which trajectories loop around many times before repeating. (D) Example of a strange (chaotic) attractor. (E) Example of a strange attractor. Unlike in A–D, the attractor here is symmetric about the line .

Fig. 4.

(A) Stratigraphy through one of the four downstream branches, with color indicating the Shields stress at the time of deposition. The solid red line indicates the vertical distance , which characterizes the predictability of stratigraphy. Parameter values are the same as those in Fig. 3. (B) Time series of bed elevation for the upstream and downstream points in the profile over the same time interval as A. (C) Time series of water-discharge fraction into the downstream branch over the same time interval as A and B. (D) Compensation statistic computed across the entire model domain. For short time windows, anticompensation occurs, whereas the system is compensational over long times (the compensation statistic decays as a power law with a slope of ). For reference, the dotted red line is the compensation time computed using a compensation length scale obtained by spatially averaging the bed-elevation SD at each point.

Dynamics

We begin by describing the types of dynamics that can be observed in our network model. At any given time, the model consists of bed elevations, where is set by the discretization of the branches. Because the system remains unchanged following an arbitrary vertical translation, in fact only bed elevation differences are needed to fully specify the system (e.g., we can subtract the upstream-most bed elevation from all others). We can then think of the time evolution of the system as a trajectory through dimension phase space. An attractor is a set of points that any point in a surrounding region of phase space (the basin of attraction) asymptotically approaches. Examples of attractors include stable fixed points, stable limit cycles, and strange (chaotic) attractors. Depending on the choice of parameters, we observe all three types of attractors in the delta-network model. It is common in other systems for the type of dynamics to depend on parameters; for instance, the Lorenz system is a well-known example of chaos but exhibits stable limit cycles or fixed points for some choices of parameters.

For the delta network, a fixed point corresponding to a symmetric partitioning at all three bifurcations and uniform aggradation across the network always exists. However, as with the single-bifurcation models discussed above, the fixed point may be stable or unstable, depending on model parameters. When the fixed point is unstable, then the flux partitioning of the network becomes dynamic. At any given moment in time, the flux partitioning tends to be uneven, causing faster sedimentation in some areas than others. However, persistent unequal sedimentation would cause the relative elevation of some areas to become increasingly low, a situation favoring increased sediment supply to those areas. Therefore, flux redistribution over time leads to uniform aggradation everywhere on average, with the rate .

For some parameter values, we observe periodic dynamics associated with stable limit cycles (Fig. 2 ). We note that these limit cycles are often asymmetric, implying the existence of another limit cycle with the asymmetry reversed (altogether, there are eight equivalent labelings of the bifurcation network; for each of these symmetries, either a limit cycle must be itself symmetric, or otherwise a reciprocal limit cycle exists). Note that even when a limit cycle is asymmetric, the average sedimentation rate over the course of a limit cycle is equal everywhere; this can be achieved through different combinations of avulsion magnitude and duration.

For other choices of parameters, we observe aperiodic dynamics (i.e., the time series of flux partitioning never repeats itself exactly). In some cases, the dynamics correspond to full avulsion, which we define as avulsion behavior where, at times, at least one branch is completely abandoned (44). Alternatively, the network may remain in the soft-avulsion regime (), in which all branches maintain a positive water discharge indefinitely. Despite the aperiodicity of the dynamics, the system remains within a bounded region of phase space, and the statistical behavior is stationary when averaged over sufficient duration. Aperiodicity does not by itself imply chaos; in fact, Fig. 2 shows an example of a quasiperiodic attractor (trajectories are confined to the surface of a 2-torus but never repeat exactly). However, in our system, quasiperiodicity appears to be rare relative to true chaos (Fig. 2 ). Next, we will show how the existence of chaos in our model can be established.

Sensitive Dependence on Initial Conditions

The hallmark signature of chaos is the sensitive dependence on initial conditions. This is defined by the exponential growth of arbitrarily small perturbations: the trajectories of two simulations with ever-so-slightly different initial conditions diverge exponentially until the difference between simulations is on the scale of the size of the attractor. In other words, given an initial separation , the distance between two trajectories is described by , where is the maximum Lyapunov exponent (MLE). The exponential divergence is not perfectly uniform across the attractor; therefore, to find , we average over many separations (57). We obtain positive values of the MLE, confirming that our system is chaotic (Fig 3 and ). In contrast, for parameter values resulting in a stable limit cycle, the MLE is zero. We observe that the MLE appears to vary smoothly across transitions from a stable limit cycle to a strange attractor () and find that at least some of these transitions follow a period-doubling route to chaos, similar to the logistic map ().

Fig. 3.

(A) Time series of water discharge through one of the four downstream branches for two independent model runs that started with nearly identical initial conditions. At first, the two time series are essentially indistinguishable, but they eventually diverge from each other completely due to sensitive dependence on initial conditions. For reference, the Lyapunov time is shown as a solid red line. (B) Distance between the two trajectories shown in A through time. Distance is measured as the L2 norm of the displacement between trajectories in phase space. The divergence between trajectories is exponential on average, with wiggles due to nonuniform divergence across the attractor. Trajectories continue to diverge until the separation is on the scale of the attractor. For reference, a straight line indicates the Lyapunov exponent of the attractor, which was obtained by averaging over many separations. (C) Plot of the correlation integral, i.e., the fraction of randomly chosen pairs of points that have a Euclidean distance of less than between them. for small , where is the correlation dimension. (D) Power spectrum obtained from a time series of discharge fraction through a single branch. For reference, vertical lines show the frequency of the unstable limit cycles obtained either when the upstream bifurcation is perfectly symmetric or when the two downstream bifurcations are perfectly symmetric.

Characterizing the Attractor

When viewed in phase space, strange attractors are typically fractal. We can quantify the fractal dimension of the attractor using the correlation dimension (23). This technique involves selecting random points on the attractor and measuring the distance between all possible pairs of points. The correlation integral is the fraction of those point pairs that have a distance between them of less than . For small , the correlation integral is a power-law function , where is the correlation dimension, a type of fractal dimension (Fig. 3). We have observed fractal dimensions in excess of four, although the fractal dimension varies as a function of parameter values. Typically, at parameter values near transitions to a limit cycle, the fractal dimension tends to be lower (slightly larger than two). The structure of these lower-dimensional chaotic attractors can be visualized using a 3D plot of the phase space, as shown in Fig. 2, whereas the higher-dimensional attractors cannot be adequately represented through three variables.

The existence of a strange attractor does not imply that avulsion is random. Although the avulsion process is not perfectly periodic, there may still be a dominant frequency. The simplest way to check for dominant frequencies is through a Fourier transform of the time series of water-discharge partitioning through one of the four branches. When we take the Fourier transform of the time series of a limit cycle, we obtain a sharp peak at the frequency associated with the periodicity of the signal, with decaying peaks at the harmonics. Compared to the limit cycle, we find that power is distributed across a broad range of frequencies in the case of a strange attractor (Fig. 3). Nevertheless, a broad peak can typically be observed in the Fourier spectrum. We observe that this peak frequency roughly corresponds to the average time between discharge peaks in the time series, but for a chaotic time series, there is enough variability in the timing between discharge peaks to make long-term prediction impossible. As a point of reference, we can find the frequency at which the upstream bifurcation avulses independently of the downstream bifurcations and the frequency at which the downstream bifurcations avulse independently of the upstream bifurcation. This is done by leaving either the upstream or downstream bifurcations unperturbed; the bifurcations are unstable in that any arbitrarily small perturbation to the symmetry grows, but without the initial perturbation, the bifurcations remain symmetric indefinitely. We do not find a consistent relationship between the dominant frequency of the strange attractor and the frequencies of the unstable limit cycles; it may line up closely with one of the unstable limit-cycle frequencies, or it may be higher or lower than either unstable limit-cycle frequency ().

Application to Stratigraphy

The lack of long observational time series presents a challenge to testing our model against field data. However, long records of delta dynamics are potentially accessible from stratigraphy. Our model does not generate stratigraphy directly, because we assume a single grain size, and we assume a fixed channel network (i.e., all avulsions are reoccupational). As a workaround, we can plot the Shields stress at the time of deposition, which should manifest itself via grain size and/or sedimentary structures. Additionally, statistics based on time variation of the bed-surface elevation are straightforward to compute using the model and have already been measured for experimental deltas and ancient systems recorded in stratigraphy (58).

We synthesize an image of a longitudinal slice of stratigraphy generated by one of the four branches, colored by Shields stress at the time of deposition (Fig. 4). Assuming that Shields stress can be considered a proxy for grain size, we observe features that are qualitatively consistent with observed stratigraphy: most beds fine upwards (59), we observe both stratigraphically transitional and abrupt avulsions (9), and some “mud plugs” (deposits formed at low Shields stress) are preserved (60).

For a quantitative comparison, we can calculate a dimensionless compensation statistic for the bifurcation network:where is the deposition rate measured at a point over the time interval , and is the long-term average sedimentation rate of the system (61). As shown in Fig. 4, the deposition rate is not constant through time and incorporates periods of erosion and/or stasis (58). However, averaged over sufficiently long time, the mean deposition rate at any given point approaches the average, . Therefore, as the measurement interval increases, decays, because more of the variability is averaged out. When the decay follows a power law, we can define the compensation index as the exponent of that power law. If , then fluctuations in the sedimentation rate are random. implies some degree of compensation; when , then the system is purely compensational, meaning that any increment of anomalously high or low sedimentation is canceled out by the next increment of sedimentation. Finally, implies anticompensation, meaning that anomalously high or low sedimentation rate is more likely than random to persist. Previous studies have shown that the compensation index is scale-dependent; at long time scales, , but lower values of are observed over short time scales (62, 63). These studies proposed that the compensation time scale determines where the transition to compensational behavior occurs, where , and is a roughness length scale of the delta surface. Consistent with previous work, we find that our model displays time scale-dependent compensation statistics. At short time scales, the system is anticompensational, and there is a gradual transition to pure compensation at long time scales. To select a relevant roughness length scale, we calculate the SD of the detrended bed-elevation time series at each point and average over the entire domain. We find that the resulting compensation time scale reasonably coincides with the gradual transition to pure compensation. An important implication of long-term purely compensational behavior is that whereas the instantaneous sedimentation rate at a point becomes unpredictable for forecast windows beyond the Lyapunov time, the sedimentation rate at a point averaged over a period longer than the compensation time scale is predictable and matches the basin-wide subsidence rate.

Discussion

Our results based on a model of a simple delta network indicate that two-way coupling between upstream and downstream bifurcations in a network can lead to chaotic dynamics in the flux partitioning of the network. This is in contrast to the behavior of a single bifurcation, which may produce limit cycles but never chaos. Although through a different mechanism, the possibility for interaction between bifurcations in a network was previously noted for braided rivers (64). In our system, the physical mechanisms behind coupling are as follows: firstly, the more intuitive mechanism is the upstream-to-downstream coupling. A change in the water partitioning in the upstream bifurcation instantaneously changes the upstream boundary conditions for the downstream bifurcations by affecting both the width-to-depth ratio and the Shields stress. As shown by previous studies (43, 51, 53), bifurcation asymmetry is sensitive to both of these parameters. Secondly, downstream-to-upstream coupling occurs because avulsions in the downstream bifurcations propagate a signal upstream through the bed-elevation profiles of their feeder channels. This alters the branch slopes at the upstream bifurcation, thereby affecting the flow partitioning. We can conceptually understand how avulsion affects the bed elevation of the feeder channel as follows: typically, following an avulsion, the new channel is steeper than the old one. This locally steeper region of the bed-elevation profile produces a higher sediment-transport rate. Via the Exner equation, this causes erosion (or at least reduced deposition) that propagates upstream. In an extreme example, following the 1855 avulsion of the Yellow River, 5 to 10 m of erosion propagated up to 100 km upstream (8). We note that backwater effects are not included in our model; therefore, the downstream-to-upstream signal propagation is purely morphodynamic, a variant on the morphodynamic backwater concept (65, 66). In contrast, the upstream-to-downstream transfer of information is effectively instantaneous, because it takes place through the water discharge. Similarly, if we were to incorporate backwater effects in our model, this would provide a mechanism for instantaneous coupling from downstream to upstream through the hydrodynamics (67).

The two-way coupling between bifurcations in a delta network is a core concept of this work. Through appropriate analysis of experiments, high-fidelity models, or satellite/historical/stratigraphic records, we believe that it should be possible to obtain evidence of this coupling in physical delta-network systems. For example, a validating experiment could be based on a fixed-width/location distributary network, ideally under flow conditions that promote bifurcation stability but limit the formation of alternate bars (48), which have been shown to complicate discharge-partitioning dynamics in single-bifurcation experiments (45). We also expect two-way coupling between bifurcations to occur under more general conditions, e.g., within delta networks where channel widths and network geometry can change through time.

Our study is not the first to use the BRT model for a system involving multiple coupled bifurcations; Kleinhans et al. (49) modeled a network consisting of 15 coupled bifurcations and confluences. In their study, the modeled river network evolved toward a frozen state with highly asymmetric flux partitionings at bifurcations. Similarly, in a Delft3D model, delta networks became frozen when they reached the edge of the computational domain, which caused deposition to cease (68). These studies highlight the importance of long-term deposition in sustaining avulsion dynamics (44). Indeed, when we run our simple delta-network model with a bypass fraction of 1, the system inevitably evolves toward a frozen state, i.e., a fixed point, with a potentially asymmetric flux partitioning.

We find that differing yet physically reasonable parameters can result in qualitatively different dynamics. Whereas the stability of the symmetric fixed point depends on model parameters such as the width-to-depth ratio in a straightforward way, the occurrence of periodic vs. chaotic dynamics is less intuitive; for instance, we find that stable periodic behavior can occur at parameter values that are book-ended by chaotic behavior. So-called “periodic windows” are common in chaotic systems. Our results suggest that a change to the forcing conditions of a delta (e.g., due to climate or tectonics) could cause a change from predictable to unpredictable behavior or vice versa.

In the case that a delta network is within a chaotic regime, what are the bounds on its predictability? The Lyapunov time () is the characteristic time scale over which a chaotic system becomes unpredictable. The Lyapunov time is not the same for all chaotic attractors produced by the model, implying that some delta networks may be predictable for longer time than others. However, we find that a Lyapunov time on the order of is typical (e.g., Fig. 2), meaning that the delta network becomes unpredictable on time scales on the order of the characteristic avulsion time. More specifically, we can ask how long it takes for the error in a prediction to grow to an unacceptable level. For example, assume that we know the initial conditions within an accuracy of , i.e., with an error of of the upstream channel depth. We can then ask at what time horizon our error grows to an unacceptable value, say . If we assume a typical Lyapunov time of 1, then the time horizon is 4.6 characteristic avulsion times. If we improve our estimate of the initial conditions by two orders of magnitude, the time horizon is improved to 9.2; the accuracy improves only logarithmically. When considering the predictability of stratigraphy, we can recast the Lyapunov time as a length scale, which characterizes the vertical distance over which information at one stratigraphic level can be used to predict the overlying stratigraphy. Because we nondimensionalize time by the characteristic avulsion time scale (channel depth divided by mean deposition rate), the time scales reported above are equivalent to length scales in terms of multiples of the upstream channel depth.

Predicting delta-network dynamics to five characteristic avulsion times or even just one is beyond the reach of existing state-of-the-art models such as Delft3D. Our simple estimate above for the time horizon at which a network becomes unpredictable suggests that there is substantial room for improvement in short to medium-term forecasting of delta networks before fundamental constraints on predictability take over. Predictions could be improved not only by increasing model resolution and better constraining model parameters but also by increasing our understanding of the relevant physical processes in order to improve the fundamental equations we choose for modeling delta dynamics.

Even when the dynamics are chaotic and therefore unpredictable in the long-term, they are far from random. Time series obtained from the model possess a strong degree of structure despite being nonperiodic. For example, we found that both unusually long and short waiting times between branch reoccupations are rare relative to random (). Given a sufficiently long averaging window, statistics such as averages, the power spectra, and range of asymmetry values can be predicted. For instance, the sedimentation rate at any point matches the basin-wide subsidence rate when averaged over a window longer than the compensation time scale. The time scale-dependent compensation statistics we obtain from the model can be contrasted with the prediction of a compensation index of if sedimentation rate were to vary randomly. We typically find that even when the dynamics are chaotic, there are typically preferred avulsion frequencies (although there is always some range; otherwise, the system would be periodic and therefore predictable). Similarly, the statistics of discharge-partitioning values can be characterized; we have found strange attractors that remain in the soft-avulsion regime (no branches are ever abandoned) and others that at any given time may have all branches open or as few as one branch open. We have even found strange attractors that are asymmetric, implying that a strange attractor with the asymmetry reversed also exists (Fig. 3). Even though chaos implies that the long-term evolution of delta networks is unpredictable, the statistical behavior is predictable (69). This is analogous to how long-term predictions of the weather are impossible, and yet one can predict the statistics of the weather (i.e., climate) on time scales well beyond the weather-predictability horizon. The compensation time scale is a measure of the time required to transition from delta “weather” to delta “climate.”

Our results show that highly complex behavior can occur even in a highly idealized system with many simplifying assumptions. In particular, we use the assumption of a fixed network geometry, implying that all avulsions are reoccupational. Although this is clearly an oversimplification of natural delta networks, a growing body of literature shows that avulsions commonly exploit previously formed channels (6, 8, 31, 70, 71). Additionally, our simulations are for a delta network consisting of just two orders of bifurcations. We have found that adding an additional order of bifurcations does not fundamentally alter the chaotic dynamics described here and produces a similar predictability horizon, although the fractal dimension is higher (). Given the simplifying assumptions of our model, and the simple delta-network geometry we chose, we expect chaos to be common in delta networks in the field, which have more degrees of freedom than our model.

Our example of a simple delta network demonstrates a particular way in which chaos can arise due to coupling between processes that are not chaotic in themselves. For instance, in the Lorenz system, it is the nonlinear coupling between the state variables that gives rise to chaos. Given that most geomorphological systems are nonlinear, dissipative, and involve many coupled processes, chaos is likely common in these systems (38). Models from previous studies may be chaotic: Murray et al. (19) found that their braided river model exhibits “apparently unpredictable changes in configuration indefinitely,” despite unchanging external forcing. In Delft3D delta simulations, varying the white-noise initial condition resulted in deltas that differ in details but whose gross-scale morphology remains similar (72). These examples are suggestive of sensitive dependence on initial conditions, even though neither formally demonstrated the existence of chaos. We propose that it is useful to know whether there are fundamental limitations to the prediction of a given geomorphic system. Whenever detailed predictions are required, we suggest that varying the initial conditions within their range of uncertainty should be included as part of a Monte Carlo analysis, in addition to the more common approach of quantifying uncertainty due to parameters. While sensitive dependence on initial conditions ultimately limits our ability make accurate predictions, estimation of the Lyapunov exponent provides a clear, quantitative measure of how we expect prediction accuracy to decay.

Conclusions

We present a simple model of a delta network consisting of an upstream bifurcation coupled to two downstream bifurcations. We find: 1) While an individual bifurcation produces periodic (nonchaotic) dynamics, two-way coupling between upstream and downstream bifurcations can produce chaos. 2) Our model produces qualitatively realistic stratigraphy and generates scale-dependent compensation statistics that are consistent with ancient and experimental delta strata. 3) The Lyapunov exponent estimated from our model implies that the loss of predictive capacity due to chaos occurs on a time scale comparable to the avulsion time. Translated to stratigraphy, this implies a loss of prediction capacity on a vertical scale of a few channel depths. 4) While chaos implies that in the long-term, delta “weather” (precise configuration) is unpredictable, delta “climate” (statistical description) can be predicted. The required time for averages to become predictable is given by the compensation time scale. 5) Based on estimates of the predictability horizon, we suggest that there is substantial room for improvement in the modeling of delta networks before fundamental limits to predictability are encountered.

Materials and Methods

Model Formulation.

We consider a simple distributary network with geometry as follows: an upstream branch, labeled , bifurcates into two smaller branches, , each of which, in turn, bifurcates into two smaller branches, where branches are downstream of branch , and branches are downstream of branch . To match typical hydraulic geometry, we choose and .

Each branch is described by a bed-elevation long profile, , where is distance ranging from 0 upstream to downstream, and is time.

First, we consider the evolution of these bed-elevation long profiles. The equation for mass conservation of sediment reads:where is the sediment flux per unit width, and is the bed porosity.

The sediment flux is calculated using a generic sediment transport formula, which reads:where is the Shields stress, is the grain size, is the submerged specific gravity of sediment, is gravitational acceleration, and are parameters, and is the critical Shields stress.

The equation for the Shields stress can be written:where is the water discharge per unit width, is the dimensionless Chezy coefficient, which we assume is constant, and is the water depth, which, assuming normal flow, is written:Thus, by substituting Eqs. – into Eq. , we can see that the evolution through time of the bed elevation profile in each branch is a function of itself and the water discharge.

Next, we seek a set of nodal conditions at each bifurcation, through which the coupling between the branches of the network occur.

Firstly, the assumption of water discharge continuity can be expressed:where is the channel width.

Next, following the approach of BRT, we allow for lateral interaction between neighboring branches occurring over a length immediately upstream of the bifurcation point, where is an order-one parameter, and corresponds to the branch immediately upstream. The elevation of these upstream cells is denoted , with corresponding to branches 2 to 7. Note that we use the prime to differentiate between the bed-elevation profile , which is a continuous function of , and the elevation of the upstream cell .

We assume that the water level at the entrance to the bifurcation is constant. This is written:where is evaluated from Eq. , in which the bed slope is evaluated at . Eqs. , , and must be solved iteratively to find the discharge entering each branch and the corresponding depths. Subsequently, Eqs. and can be evaluated to find , i.e., the sediment flux leaving the upstream cell and entering the branch at .

Upstream cells are allowed to exchange sediment with their neighbor. This transverse sediment flux is computed:where from continuity, , and is a parameter between 0.3 and 1.0 controlling how strongly the lateral bed slope influences the transverse sediment flux.

The bed elevation of the upstream cells evolve through time according to mass conservation:

Note that in the above equations, is the downstream boundary to branches i = 1 to 3. We assume that the slope is calculated using the average of and , allowing the sediment flux to be computed from Eqs. –.

Next, we consider the downstream boundary to branches i = 4 to 7. Following the approach of SPV (44), we introduce the bypass fraction . The bypass fraction is the fraction of the input sediment flux that exits the downstream boundaries in the long-term average. Note that this balance does not necessarily occur in the short term. At the downstream boundaries, we prescribe a Robin boundary condition on the sediment flux, written:where is defined as the topset area:For illustration, we examine the behavior of () for two end members. Firstly, in the end member, we see that the equation reduces to at , implying a deposition rate of 0. In contrast, if , the equation reduces to at , i.e., all of the sediment is trapped within the aggrading delta topset.

Finally, the upstream boundary to the entire system can be specified through the water discharge, sediment flux, and grain size, or, equivalently, the width-to-depth ratio, Shields stress, and slope.