Entropy and optimality in river deltas.

The form and function of river deltas is intricately linked to the evolving structure of their channel networks, which controls how effectively deltas are nourished with sediments and nutrients. Understanding the coevolution of deltaic channels and their flux organization is crucial for guiding maintenance strategies of these highly stressed systems from a range of anthropogenic activities. To date, however, a unified theory explaining how deltas self-organize to distribute water and sediment up to the shoreline remains elusive. Here, we provide evidence for an optimality principle underlying the self-organized partition of fluxes in delta channel networks. By introducing a suitable nonlocal entropy rate ([Formula: see text]) and by analyzing field and simulated deltas, we suggest that delta networks achieve configurations that maximize the diversity of water and sediment flux delivery to the shoreline. We thus suggest that prograding deltas attain dynamically accessible optima of flux distributions on their channel network topologies, thus effectively decoupling evolutionary time scales of geomorphology and hydrology. When interpreted in terms of delta resilience, high nER configurations reflect an increased ability to withstand perturbations. However, the distributive mechanism responsible for both diversifying flux delivery to the shoreline and dampening possible perturbations might lead to catastrophic events when those perturbations exceed certain intensity thresholds.

River deltas are depositional landforms forming downstream of major rivers, often home to large populations and important natural resources (1–10). In the last decades, many deltas of the world have been under threat from a range of stressors, including sea-level rise, upstream dam development, and local exploration (2, 4, 5, 11–17). Deltas are nourished by channel networks whose connectivity constrains, if not drives, the evolution, functionality, and resilience of these systems. Remarkably, the properties of delta channel networks differ substantially from the tree-like topology of the rivers that feed them (18). Tree-like networks, defined by the absence of loops, are characteristic of tributary river networks and are found abundantly in nature across different systems and scales (e.g., botanical trees, veins of leaves, blood vessels, lightning, and river networks). The propensity of nature in choosing tree-like configurations has been grounded as an optimality principle. Specifically, tributary river channel networks achieve minimal total energy dissipation, that is, minimal loss of potential energy as water and sediment flow downstream, albeit often manifesting as feasible optimality, that is, a dynamically accessible local minimum due to initial conditions and other constraints (19–22). Similar to river networks, vascular networks in biological systems (e.g., animals, plants, insects, etc.), which transport materials through space filling fractal networks of branching tubes, achieve states of minimal energy dissipation (23, 24). Analogous optimality principles have also been suggested to constrain processes as diverse as root water uptake in plants (minimization of internal dissipation) (25) and land surface energy and water balance (maximization of power) (26–28).

Despite the importance of deltaic systems and the recent advances in quantifying their connectivity properties (18, 29–34), an optimality principle for the organization of their distributary channel networks, akin to that existing for the tributary networks, remains elusive. A recent framework based on spectral graph theory (29, 30, 32) sheds light on the topologic and steady-state flux partitioning characteristics of delta channel networks and their relationship to underlying morphodynamic controls (e.g., sediment composition and tidal and wave energy), paving the way for quantitative delta classification and inference of process from form. However, the diversity of topologic structures of channel networks across a broad spectrum of deltas, and the fact that deltas are highly dynamic systems within which topology and flux partition coevolve, make it challenging to find a universal first-order optimality principle (e.g., minimization of energy, maximization of entropy, or minimization of free energy) governing their formation.

We are physically motivated by the foundational principle that deltas build land by spreading their fluxes on their delta top, as opposed to creating single pathways to the ocean which would diminish the formation of islands that retain sediment and nutrients and reduce land-building potential. This notion resonates with previous results applied to tidal deltas showing that tidal channels self-organize to uniformly distribute the tidal prism across the delta (35). Under this premise, we postulate the existence of an optimality principle in delta channel networks in terms of achieving configurations that maximize the diversity of flux delivery from any point of the network to the shoreline. The time scale of topologic reorganization in deltas, which is mainly through major channel avulsions, is of the order of decades to millennia, depending on the deltaic system (36, 37), whereas the flux adaptation on a given channel topology has a characteristic time scale which is orders of magnitude smaller. This separation of scales allows us to study the system on a fixed topology and focus on whether an optimality principle governs the distribution of fluxes on that topology.

Given the nature of our hypothesis, that is, maximization of the diversity (or uncertainty) of flux delivery from any point of the network to the shoreline, we adopt an information theoretic approach to quantify uncertainty based on Shannon entropy (38). Note that Jaynes (39) introduced a formalism demonstrating the equivalence (under thermodynamic equilibrium) of the statistical mechanics and information theory approach to entropy. Furthermore, it has been argued (not exempt from certain controversy) that the mathematical framework formulated by Jaynes (39) serves as a generalization of the statistical mechanics framework for both equilibrium and nonequilibrium systems (40, 41). In this work, we propose the notion of nonlocal entropy rate () and suggest by comparative analysis of field and numerical deltas that indeed deltas self-adjust their flux partition so that they maximize their flux diversity. Note that by “maximization” we do not refer to a global optimization, that is, deltas generally do not achieve flux configurations that correspond to the absolute maximum value of but exhibit configurations corresponding to dynamically accessible (local) maximum (feasible optimality). We also discuss the possible implications of the proposed optimality principle with respect to delta resilience in response to perturbations, arguing that flux distributions characterized by extreme values of are more resilient, in that a local perturbation (e.g., flux reduction in a channel) will affect the least the distribution of fluxes at the shoreline outlets.

Nonlocal Entropy Rate (nER)

Entropy quantifies the uncertainty in the occurrence of events (38), that is, the amount of information needed to describe the outcome of an experiment. Uncertainty intuitively emerges from the notions of probability and surprise. For instance, given a discrete stochastic process , such as rolling a six-sided die, with specified probability distribution of outcomes, for example , the occurrence of rolling a 3 when all sides are numbered 3 produces zero surprise. Conversely, on the same die rolling a number other than 3 would produce infinite surprise. Mathematically, this surprise is defined as . Therefore, the uncertainty of an event, , is the product of times its probability of occurrence . Thus, either the occurrence of a completely certain event ( = 1) or an impossible one ( = 0) introduces zero uncertainty ( = 0, by convention  = 0). The uncertainty is maximal when all of the possible outcomes have the same probability of occurrence 1/ (e.g., fair die where the number on each face has a probability of occurrence 1/6). The total uncertainty or entropy, , of a set of discrete outcomes with probabilities is equal to the sum of the uncertainties corresponding to each outcome (42):We aim to develop an entropic metric for delta channel networks that quantifies the diversity of flux delivery to the shoreline. Thus, a delta channel network with low entropy, that is, with low uncertainty in flux pathways, would be one with a dominant channel that carries most of the flux. This configuration is inherently unstable, assuming no redistribution of sediment by marine processes, because that channel would prograde until its water slope reduces and the channel avulses down a steeper path (43, 44). However, a delta channel network that achieves a configuration where all of the possible pathways that drain water and sediment from any point of the network to the shoreline are equally probable would have the highest entropy. This uncertainty in delivery would have the stabilizing effect of spreading sediment evenly across the shoreline. Notice that maximizing the uncertainty of water and sediment flux delivery to the shoreline does not necessarily imply networks that proportion fluxes equally at every single bifurcation. In fact, this would be inconsistent with morphodynamic theories (45, 46) which showed that asymmetric flux partition at a junction is a requirement for stability.

Conceptualizing delta channel networks as graphs, where nodes correspond to junctions or bifurcations, and links represent channels (), we define a metric of uncertainty of water and sediment flux pathways from the node to the shoreline aswhere is the transition probability from node to outlet node . Alternatively, represents the fraction of water and sediment flux from node that eventually drains to outlet (see for details on the computation of ). We refer to as the nonlocal entropy for node , to emphasize that the transition probabilities are between nodes and the shoreline nodes as opposed to among neighboring nodes (local). This notion acknowledges important nonlocal effects on delta dynamics, such as the hydrodynamic backwater where the water surface slope is dependent on the water depth at the shoreline and the local slope between subsequent bifurcations (36, 44).

By weighing with the normalized steady-state flux at each node , we define the of the complete delta channel network:where is a normalization factor computed for each node and represents the number of outlet nodes that can be reached from node . The normalized nER admits values in the interval [0,1].

Results and Discussion

We hypothesize that deltas distribute the flux at each bifurcation to maximize the uncertainty in the delivery of fluxes from any point of the delta to the shoreline, that is, to achieve a dynamically accessible maximum of . We computed the for 10 deltas with diverse morphodynamic environments and channel complexity (Fig. 1; see for the extracted channel networks and physical information about the deltas). We compared the computed for each field delta (based on the actual flow partitions) against randomizations of the flux partition at each bifurcation [sampled from a uniform distribution in the interval ] holding the network structure constant. Despite the broad range of climate, discharge, and sediment influencing these deltas, all but one (the Niger delta) have flux configurations that exhibit extreme values of , that is, 9 of 10 field deltas have above the 90th percentile of the random distribution (probability of exceedance < 0.1) (Fig. 2). The fact that the computed value of does not correspond to the absolute maximum of the distribution is not surprising. Natural systems have been argued to achieve stationary configurations that do not correspond to the absolute optimum of the functional describing their organizational principle but to local optima that are accessible given the initial conditions, constraints, and the system dynamics, known as the feasible optimality principle (19–22). Beyond field deltas, we also applied the analysis to numerically simulated deltas that formed under varying incoming sediment grain sizes (32, 47, 48) using the physically based hydromorphodynamic model Delft3D ( for further details). The results show that five of the six numerical deltas exhibit extreme values of with probability of exceedance  < 0.1 (Fig. 3), further supporting our optimality hypothesis. Note that the simulated delta ( = 0.01 mm) that does not satisfy the optimality of is an extreme case, in terms of cohesiveness, for field deltas. Very cohesive banks are harder to erode and form levee breaches infrequently, delaying the triggering mechanism of avulsions. As a result, the system maintains itself at states at which the fluxes are not at equilibrium with its underlying channel network topology. This is reflected in suboptimal states of flux distribution and thus .

Fig. 1.

Field deltas and their corresponding location. Ten deltas with diverse morphodynamic environments and of various degrees of channel complexity were analyzed in terms of their nER (see for further details about the deltas). Satellite images provided by Landsat/Copernicus, NASA, Digital Globe, and CNES/Airbus were extracted from Google Earth. We acknowledge their respective copyrights.

Fig. 2.

for 10 field deltas. Green stars represent the values of nER computed for each field delta, using channel width (extracted from Landsat) as proxy for flux partition in bifurcations. We compared the values of for each delta with 105 randomizations of flux partitions (histograms). Nine out of the 10 deltas analyzed exhibit a maximal value of , defining maximal as a value where the probability of exceedance, , by a random realization is less than 0.1.

Fig. 3.

for simulated deltas. We examine the nER of numerically simulated deltas obtained by the Delft3D model. The simulated deltas are river-dominated, with no vegetation, and with a lognormal distribution of incoming sediment size with median grain size varying from 0.01 to 1.0 mm and the same variance in the log space (for more details see ). These deltas exhibit a wide variety of channel network topologies as shown in ref. 32. Similar to the analysis conducted for the field deltas, green stars represent the values of nER using channel width (extracted from simulations) as proxy for flux partition in each bifurcation. Compared with randomizations of flux partition (histogram), five out of six deltas analyzed exhibit a maximal value of , defining maximal as a value where the probability of exceedance by a random realization is less than 0.1.

Our results suggest that the flux partitions at each bifurcation, which have evolved naturally, are not random but rather follow a rule that optimizes the delta system as a whole. In fact, an interesting paradox arises from our analysis. Although the entropy introduced locally by each bifurcation, considered as an insulated unit, is suboptimal [the maximum would correspond to a symmetric bifurcation which is not consistent with stability theory of delta bifurcations that requires asymmetric local flux partition (45)], the specific assemblage of those bifurcations forming the delta network as a system is optimal (in terms of ) and consistent with maximization of the diversity of fluxes delivered to the shoreline.

Turning attention to delta dynamics, we further hypothesize that during an avulsion the delta would decrease (see schematic in Fig. 4). This is because during this phase of topologic reorganization the flux distribution inherited from the previous channel network structure is in general suboptimal with respect to the incipient channel network reworked by the avulsion (i.e., the new channel structure created during the avulsion received a disproportionally small share of the flux, creating an asymmetry in the flux delivery to the shoreline and thus reducing the value of ). Testing this hypothesis in field deltas is challenging because avulsions occur infrequently. However, using numerical models we can observe that during an avulsion cycle the drops significantly at the onset of a new flow path and following the abandonment of the old channel (Fig. 4). Since the time scale of the avulsion itself is negligible in comparison with the lifespan of the topologies before and following the avulsion, it is observed that the flux partition is able to self-organize to achieve a configuration that maximizes . This supports our assumption that the time scale of the flux reorganization is several orders of magnitude smaller than the characteristic time scale of topologic reorganization which is set by the time lapse between avulsion cycles.

Fig. 4.

nER during an avulsion cycle. (A) The results presented in Figs. 2 and 3 support the hypothesis that river deltas self-organize to maximize nER. However, deltas in dynamic equilibrium experience avulsions, periodically causing major reorganization of their channel networks. Our hypothesis is that during these periods of major reorganization delta nER drops significantly until the delta fluxes are self-organized back to a new structure as shown in the schematic of A. (Insets, from left to right) Schematics of the five stages within a single avulsion cycle (notice that for illustration purposes, only the main channel was drawn to better depict the avulsion cycle): (I) starting from a given channel network topology, (II) a new channel structure is created downstream of the avulsion node draining initially a small part of the flux, but progressively receiving larger proportions of fluxes until (III) the new and the old channel structures receive similar amounts of flux, and eventually transition to (IV) a configuration wherein the new channel structure has more flux, to finally lead to (V) the abandonment of the former channel structure downstream of the avulsion node. (B) We have tested this hypothesis using Delft3D simulations and analyzing subsequent instances of a delta evolution, wherein an avulsion cycle is observed. The probability (Prob.) of exceedance of the delta nER by a randomization of the flux configuration is displayed. The results agree with the posed hypothesis, showing that nER drops at the beginning and end of the avulsion cycle, keeping a maximal value otherwise.

An important implication of this optimality principle can be interpreted in terms of the resilience of deltas to withstand perturbations. Intuitively, if a perturbation (e.g., flux reduction) is applied to a delta during its high- state, the perturbation will be damped as it will spread through the diverse pathways connecting the delta top to the shoreline. However, if the same perturbation is applied to a delta in a low- state, the perturbation will be more confined to a localized part of the delta but will exert a more severe disturbance. As revealed by our analysis, river deltas operate in configurations characterized by high values of , supporting the idea that deltas self-organize to achieve resilient morphologies priming self-maintenance. As a word of caution, especially relevant in the current scenario where deltas are subjected to increasing anthropogenic stresses, this distributive mechanism that dampens the intensity of perturbations can also lead to delta-wide catastrophic disturbances and tipping points when those perturbations exceed certain thresholds.

Conclusions

Deltas are highly productive regions supporting extensive agriculture and aquaculture and diverse ecosystems and containing natural resources such as hydrocarbon deposits. Climate change and human actions, both in the upstream basins and locally, act as stressors on these landscapes, calling for a thorough understanding of these complex systems and their response to perturbations. We examined the existence of an optimality principle that governs the self-organization of water and sediment fluxes on delta channel networks. Specifically, (i) we put forth the hypothesis that maximizing nER, which quantifies the diversity in flux delivery to the shoreline, is a selective criterion in the evolutionary dynamics of delta networks; (ii) we tested this hypothesis by analyzing 10 field deltas of diverse complexity, age, and environmental settings and showed that all but one, the Niger delta, exhibited maximum nER; (iii) we further supported the existence of an optimality principle by analysis of Delft3D simulated deltas; (iv) we showed that during major reorganization, such as avulsions, nER exhibits suboptimal values and increases back to high rates (maximum values) when the flux distribution self-adjusts to the new delta channel network topology; and, finally, (v) we discussed the relation between entropy and resilience, arguing that delta flux configurations characterized by maximal nER are more resilient in the face of random perturbations. In the anthropocene where human activities have become a major agent of geomorphic change, understanding delta self-organization within an optimality perspective offers new ways of thinking about delta dynamics and disturbances that might hinder self-maintenance.

Materials and Methods

Deltas as Directed Graphs.

Tejedor et al. (29, 30) presented a rigorous framework based on graph theory within which a delta channel network is represented by a directed graph, that is, a collection of vertices (bifurcations and junctions) and directed edges (channels in-between vertices, where the direction is given by the flow). All information about network connectivity and directionality of the flow can be stored in a sparse matrix called the adjacency matrix, . Specifically, is an matrix, where is the number of vertices, and whose entry is unity if vertex receives fluxes directly from vertex (i.e., vertices and are connected by a link directed from to ) and zero otherwise. From we can derive an important matrix called Laplacian, which is equivalent to a diffusivity operator in a graph. To construct the Laplacian we need first to introduce the degree matrices for directed graphs. The in-degree (out-degree) matrix () is an diagonal matrix whose entries depict the number of links entering (exiting) vertex and are computed as the sum of the entries in the -th row (column) of . The Laplacian matrix, (), is defined as (). Tejedor et al. (29, 30) showed that certain eigenvectors of the Laplacian operator contain important topologic information of the deltaic network. Furthermore, information about flux propagation can be obtained if is replaced with the weighted adjacency matrix , where the weights correspond to the fraction of flux in link () with respect to the flux in vertex . Similar to the at-station hydraulic geometry relationship (49)—width to landscape-forming discharge—reported for tributary rivers (50) and tidal channels (51), we assume the flux partition at the bifurcation to be proportional to the width of the downstream channels (18). Note that even though we do not consider explicitly in the computation of steady-state fluxes relevant processes such as water–sediment interchange between channel and islands (31, 52, 53), vegetation (54), tides (35), and so on, all of these processes set the hydrogeomorphic attributes of the channel network. Therefore, the computation of steady-state fluxes in the channel network based on physical attributes such as channel widths can be interpreted to a certain degree as the result of the integrative effect of all of the main processes acting on a delta. For the purposes of this paper, there are two probability distributions that can be computed by simple algebraic manipulation in the above-mentioned operators, namely, the steady-state flux distribution and the node-to-outlet transition probability distribution.

Steady-State Flux Distribution.

Having a delta represented as a directed acyclic graph (DAG) allows us to compute the steady-state flux by assuming conservation of mass. For instance, Tejedor et al. (29) showed how the steady-state flux can be formulated as an eigenvalue–eigenvector problem of the Laplacian matrix of the graph. Here, we present the input/output model as a more intuitive way to compute the delta steady-state flux. Let us consider a DAG fed from the most upstream node (apex) with a constant unit flux. Without loss of generality, and for simplicity in the subsequent derivations, the apex is assumed to be labeled as node 1. Then, we can define the stationary distribution, , aswhere the column vector corresponds to the initial state, and corresponds to the weighted adjacency matrix of the graph. Alternatively, represents a constant inflow at the apex normalized to 1, and the resulting downstream response at time/distance specified by . For a DAG, there exists at least one indexing of the graph such that each offspring vertex has a higher index than its parent vertex. In this indexing, by construction, the matrix is strictly upper triangular, and therefore the matrix is nilpotent, guaranteeing the convergence of the sum as expressed in Eq. . In we prove that the stationary flux distribution can also be obtained as the stationary probability of a Markov process .

Node–Outlet Transition Probability Distribution.

The node-to-outlet transition probability, , is defined as the probability that a package of flux at node drains to outlet . Thus, the transport from each node , to the different outlets , can be understood as a discrete stochastic process with probability distribution , with . Tejedor et al. (29) showed that when a delta channel network is represented by a directed acyclic graph with a weighted adjacency matrix ,

The null space of the weighted in-degree Laplacian for the reverse graph has dimension (multiplicity of the eigenvalue zero) equal to the number of outlets ;

There exists a unique basis , , of this null space in (i.e., the basis consists of vectors each having entries) with the propertyThat is, the entry of the vector corresponding to outlet is one, and zero at all other outlets.

The value represents the portion of flux at the vertex that drains to the outlet , that is, .

Thus, if we define a matrix , whose columns form the basis of the null space of , , then the -th row corresponds to the probability distribution .