Competing feedback in an idealized tide-influenced delta network

The morphodynamic evolution of river deltas is intimately tied to flow and sediment partitioning at bifurcations. In this work, the long-term equilibrium configuration of a simple delta network is investigated by means of an analytical model, which accounts for the effect of small tidal oscillations. Differently from individual bifurcations where tidal action is always a stabilizing factor, in the case of a tree-like delta with multiple bifurcations a dual response emerges. Specifically, depending on the values of four reference parameters functions of tidal amplitude, upstream flow conditions, and channels geometry, tides can either promote or discourage an unbalanced discharge distribution. This behavior primarily concerns the apex bifurcation, which is affected by the variations of the relative tidal amplitude at the internal nodes. In turn these variations depend on how flow and sediment are diverted upstream. The stability of steady-state solutions is found to be governed by the sign and magnitude of the slope asymmetry between channels. This work provides a basic modeling framework for the interpretation of the autogenic response of multiple coupled micro-tidal bifurcations, which can potentially be extended to include in a unified scheme erosional and depositional processes typical of fluvio-deltaic systems.


Introduction
We are concerned with the morphodynamics of channel bifurcations, a representative geomorphological feature controlling the downstream flow routing in a variety of environmental settings such as coastal plains, mountain valleys, coastal wetlands, submarine fans, and lava flows [1,3,19,38,43,74]. Over the last three decades, river bifurcations in fluvial environments have been investigated through a breadth of different approaches, including theoretical [8,59,73,83], experimental [6,7,79], numerical [24,41,72], and field investigations [5,12,44]. These studies showed the natural tendency of bifurcations to present an 1 3 unbalanced flow distribution, which is reflected by a marked asymmetry of channel widths and bed elevations in the downstream anabranches [84].
In general, such asymmetric behavior can be driven by different processes. "Allogenic" factors, as in presence of an upstream evolving meander bend or migrating bars, may force the bifurcation towards an unbalanced state [60]. However, the unbalanced distribution may also arise as a consequence of an "autogenic" instability mechanism [59], whereby for particularly wide and shallow channels, the bifurcation triggers the upstream formation of damped bed oscillations taking the form of alternate bars, which steer water and sediment fluxes towards the most-carrying branch.
In order to reproduce in a simplified way the complex flow-routing processes that characterize real world bifurcations, morphodynamic modeling involved mainly one-dimensional schemes with a main channel dividing into two smaller downstream branches. In these models, the adoption of a suitable law for the partition of water and sediment fluxes is required. The two-cell model proposed by Bolla Pittaluga et al. [9] (hereafter, BRT) showed to successfully reproduce the basic elements acting in bifurcation dynamics, as demonstrated by numerical and laboratory experiments [6,24,72]. In particular, the BRT model showed that in the absence of (external) allogenic factors, the morphodynamic evolution of a geometrically symmetric ("free", sensu [59]) bifurcation is governed by a single parameter, R [8]: where * 0 and 0 are the Shields number and half width-to-depth ratio of the main feeder channel. The parameter is an order-one parameter physically representing the extension of a transverse bed elevation gradient induced in the feeder channel by the bifurcation [59]. The sediment flux is deflected with respect to the water flow direction due to a gravitational effect, whose magnitude is controlled by the so-called Ikeda parameter r [34,52].
However, in multi-thread environments like deltas and braided rivers bifurcations are not isolated elements but rather part of an interactive system [53,70,78]. Taking the viewpoint of complex networks [45,75], competing interactions break-up the possibility of understanding the whole behavior of the system just by a superimposition of the single bifurcations responses. The question then arises of what kind of physical processes are responsible for such interactions. Clearly, the flux partitioning in the upstream bifurcations affects each subsequent downstream bifurcation by changing flow and sediment transport conditions in the connecting channels. Nonetheless, this upstream-to-downstream effect is not sufficient to generate a coupling if a downstream-to-upstream response is absent. In such conditions, where the downstream branches are passively subjected to the incoming fluxes, the bifurcation is said to be "upstream controlled" [56]. Differently, if the response of the downstream reaches depends on the evolution of the bifurcation, which in turn is affected by variations of downstream conditions, a two-way coupling is established. For example, in low relief systems such as deltas, depositional landforms arising from rivers that end at the shoreline often originating a network of channels ( Fig. 1), downstream effects like a nonuniform sedimentation [63,64] or marine processes like tides [13,31,35,57,66] influence the response of bifurcations through variations of the channels slope.
In a series of recent works, Salter et al. [63][64][65] showed that a net-depositional system is able to generate a switching in the discharge partitioning over time. The observed oscillations are purely autogenic, inasmuch they derive from an internal feedback mechanism: if the fluxes are distributed unevenly by a bifurcation, the slope of the channel receiving , more sediment is found to drop through time with respect to the other channel till eventually causing the switching [21]. Interestingly, a (deterministic) chaotic dynamics is shown when multiple bifurcations are coupled in a simple delta network [65]. A similar behavior, in which feedback mechanisms are internally generated has been observed in a bifurcationconfluence system (i.e. a loop), where a river splits into two anabranches rejoining downstream [56]. In this case, the feedback is generated by a differential water surface elevation in the downstream branches, which is a function of the amount of water and sediment delivered at the bifurcation, resulting in an asymmetric energy slope affecting discharge partition upstream.
In this paper, we aim at understanding how tides potentially affect the long-term equilibrium state of river deltas. To achieve this goal, we pursue the one-dimensional modeling of an idealized delta network made up of multiple bifurcations with branching channels flowing into a tidal sea. Tidal influence is considered small enough so that the morphodynamics of the delta is primarily controlled by fluvial processes. We explicitly exclude wave activity that, basically through an alongshore sediment transport [47], is responsible for several morphological implications including the suppression of bifurcations near the coastline [36,49]. The analysis is motivated by recent theoretical advances on individual tide-influenced bifurcations [57] showing that tidal action causes a morphodynamic response akin to the one observed by Salter et al. [63] in a net-depositional system: tide deepens downstream channels through erosion, leading the branch carrying the smaller amount of water and sediment to increase its slope over time with respect to the other branch. However, this behavior might not be univocal when bifurcations are combined in series in a tree-like structure inasmuch the response of upstream bifurcations is dependent from water surface variations at the internal nodes of the network that, in turn, are affected by flow and sediment partitioning upstream.

Stability of individual bifurcations in micro-tidal environments
In a recent paper, Ragno et al. [57] (hereafter, RTB) included in the BRT model the effect of (monochromatic) tidal oscillations. Under the assumption of tidal amplitudes a * (from hereinafter the superscript denotes a dimensional quantity) much smaller than the mean flow depth in the upstream channel D * 0 , RTB investigated the (tidally-averaged) equilibrium configurations and stability conditions of an idealized bifurcation building on the previous work of Seminara et al. [69] for single river-dominated estuaries. In Fig. 2a, an example of equilibrium diagram is shown in terms of the discharge asymmetry index Q [7]: where the suffixes 1,2 and 0 denote the two downstream bifurcates and the upstream channel, respectively. In the absence of tides (black lines in Fig. 2), the model outcome is essentially the same obtained by BRT. The bifurcation parameter R is found to control the number of equilibrium solutions and the stability of those equilibria: if R is higher than a critical threshold R C , which depends on the friction and the transport formula adopted [60], flow and sediment fluxes are equally distributed at the bifurcation; on the contrary, when R is lower than the critical threshold the balanced solution becomes unstable and the system attains a stable state where one of the two branches captures most of discharge where, at the same time, the gravitational pull deviates the sediment flux to prevent the closure of the penalized branch (i.e. carrying the lowest discharge). Fig. 2 a Equilibrium diagram of micro-tidal bifurcation according to the RTB model. The discharge asymmetry index Q (Eq. 2) is plotted against the bifurcation parameter R (Eq. 1), for different values of the scaled tidal amplitude . With R ∞ C it is denoted the threshold value as obtained by [8] in the purely riverine context. Lower panels show results from the linear stability analysis, where critical threshold R C is shown for increasing values of the scaled tidal amplitude , which spans from = 0 (i.e. absence of tide, black dashed line) to = 0.8 , as a function of b 0 ( L = 1 ) and c L ( 0 = 10 ). The value of R separating the stable/unstable solution decreases as long as the tides is increased for sufficiently high/low values of 0 /L When tides are taken into account, the critical value R C at which the balanced solution becomes unstable gets smaller (Fig. 2a), thus indicating the stabilizing role conducted by tidal action. When R exceeds the threshold R C (i.e. the main channel is particularly wide and shallow), a greater fraction of water and sediment fluxes is steered towards one of the two branches. In this case the tidal action tends to re-equilibrate the bifurcation by increasing the slope of the penalized branch, therefore reducing discharge inequalities.
Through a linear stability analysis, RTB have shown that in addition to the bifurcation parameter R (Eq. 1) three other parameters turn out to play a role, namely the scaled tidal amplitude , a parameter measuring the relative importance of flow inertia 0 , and the (dimensionless) length of the downstream branches L: where * is the frequency of the tidal wave, L * the (dimensional) length of the branches, L * B0 = D * 0 ∕S 0 is the so-called backwater length [50], U * 0 , D * 0 and S 0 are reference values for velocity, depth and slope corresponding to normal flow conditions, which are assumed to be met in the upstream channel.
For a given value of , the magnitude of the stabilizing effect increases as the branches are shortened and/or 0 is amplified (Fig. 2b-c). Interestingly, if the tidal influence is sufficiently strong there is the possibility for a configuration in which discharge is equally divided even in the limiting case R → 0 . Such conditions are expected to be met in fine-grained sand-bed rivers, which are typically characterized by high values of the Shields number. In these systems, most of sediment is carried in suspension and therefore is weakly deflected by the gravitational pull at the bifurcation [37,60]. Thus, according to Eq. (7) the sediment flux divides proportionally to the water discharge and, in the case of a purely fluvial (i.e. absence of tides) bifurcation, this partitioning would necessarily lead to an unbalanced configuration. However, for micro-tidal bifurcation the theoretical model suggests that for a given value of the parameters 0 and L there exists a "tidal threshold" that allows for a stable balanced bifurcation. In this case, the stability of the system is purely "downstream-controlled" by tidal action.

A simple model for the equilibrium of tide-influenced delta networks
In Sect. 2 we have seen that in the case of a single bifurcation tides exert a stabilizing effect that acts to prevent the closure of one of the two branches in the long-term. However, deltaic networks are formed by several interconnected bifurcations that can interact with each other [20,25,29,65,77]. Let us consider an ideal delta network, which is assumed to be composed by k orders of bifurcations (sensu [23]). As a consequence, the number of nodes (N) of the system is equal to: Following RTB, each branch of the network is assumed to be a straight estuary, where the banks are fixed (i.e. the planform geometry does not change in time) and the bottom is composed by a uniform sand with a reference grain size d * . Seaward, debouching channels are subject to monochromatic tidal oscillations with a constant amplitude in the assumption ≪ 1 . Landward, the main upstream channel is fed by a constant discharge Q * 0 , with the flow carrying a constant sediment flux Qs * 0 in equilibrium with fluid flux. The geometry of the network (i.e. lengths * and widths W * ) is assumed to be symmetric about each node. Notice that in the present problem there are no forcing effects potentially leading to a break-out of the symmetry. In other words, any externally-imposed slope advantage that may alter downstream rating curves is not considered.
Channel widths are assigned following an hydraulic geometry law empirically derived from delta networks. Specifically, we rely on the relationship derived by Andrén [2], which has been employed and tested against field and modeling studies [22,23]: where W * 0 is the width of the upstream feeder channel. Analogously, each subsequent bifurcation is placed at a distance that follows a power-law distribution [22,23,71]: where L * is the length of the first-order branches (i.e. k = 1 ). We note that the choice of different geometric laws for widths and lengths does not change the qualitative behavior of the model.
For modeling the response of the bifurcation we employ the previously introduced two-cell BRT model. The nodal relationship proposed by BRT for the transversal exchange of water ( Q * y ) and sediment ( Qs * y ) for a generic node of order k, taking into account (5), reads as follows: with * the bed elevations of the two cells, W * 0 long, which are able to exchange water and sediment upstream the bifurcation (for a detail of the two cells, see Fig. 3). Following RTB the bed elevations in the two bifurcates are affected by tidal motion, which induces a correction with respect to the normal flow of O( 2 ) . This correction turns out to be a function of the relative tidal amplitude at the downstream nodes (or the sea mouth in the case of the terminal branches), channels length, channels slope, and on the flow depth values in normal flow conditions. Taking advantage of the fact that is a small parameter, it can be shown that the water surface elevation and flow depth at each node have the following expressions [69] (see for details of the derivation [10,57]): , D u is the dimensionless uniform flow depth, and the expression for the coefficients Ψ , Θ , and are given in the Appendix 1.
In order to determine flow and bed topography in each channels, five internal (i.e. matching) conditions at the bifurcations node are required: the conservation of water and sediment mass, the equality of the water surface elevation for each cell, and the BRT relationship (7). At the end, for each node the following algebraic nonlinear system of equations is obtained: where the variables are made non-dimensional in the form: 3 Plan view of the idealized micro-tidal network composed by k = 2 orders of bifurcations (sensu [23]) hence consisting of a total of seven branches coupled through three bifurcations. The geometry of the network is symmetric at each node, with channel lengths and widths assigned following the empiricallyderived relations (5) and (6) It is worth noting that tidal amplitude in the internal nodes (see Eqs. 8a-b) is not assigned a priori, but depends on the amount of water and sediment delivered by the channels, and it can be computed through the following expression (see Appendix 1 for the details about the derivation): The inequality in (11) is a direct consequence of the problem structure: if k = 1 , the single bifurcation case analyzed by RTB is obtained. Moreover, for the terminal branches of the network the scaled tidal amplitude coincides with the one assigned at the seaward boundary ( ). As expected, | N z−1 tends to vanish over distances that are longer compared to the backwater length of the single branches. Finally, it is important to underline the steady character of the present analysis. This means that if each branch of the network is supposed to be in equilibrium (sensu [69]), this is equivalent to assume a much slower timescale for the evolution of the channels with respect to the timescale of the bifurcation evolution, which in turn is slower than the typical tidal period. This hierarchy of timescales allows for studying the response of the network as a sequence of quasi-equilibrium states, where tidal fluctuations are accounted for just by changing the amplitude .

Contrasting response of a delta network to tidal action
In analogy with the work of [65] and from observations showing that tide-influenced deltas are often composed of a limited number of distributaries [49,76], we consider a simple network made up of two orders of bifurcations (i.e. k = 2 → N = 3 , see Eq. (4)). We have noted that an additional order of bifurcations does not change the qualitative response of the system and the key drivers governing its stability. Thus, at the end the problem requires the solution of a nonlinear system of twelve algebraic equations, which can be cast in terms of the channel slopes and (dimensionless) normal flow depths of each channel, once suitable closure formulas for the sediment transport and friction are employed. Sediment fluxes are computed by means of the Engelund & Hansen [26] formula: where is the relative submerged sediment density and C f is the friction coefficient, here expressed through the Manning-Strickler formula under the hypothesis of infinitely wide channels: with n * the Manning coefficient. Consequently, the dimensionless water discharge is found to be proportional to Q ∝ D 5∕3 u √ S∕S 0 . The model requires as input just few parameters, namely 0 and R for the reference flow in the upstream channel, the scaled tidal amplitude at seaward boundary , and the dimensionless length L of the first-order branches.
Firstly, the attention is devoted to the behavior of the apex node N 0 . The equilibrium diagram (Fig. 4a) shows a striking novel feature with respect to the single bifurcation case: the response of the system to tidal action is not univocal. Indeed, the critical threshold at which the unbalanced solutions appear can either increase or decrease for different values of the tidal amplitude depending on the relative length of the channels L. This twofold behavior can be explained if results are analyzed in terms of the slope asymmetry index S (Fig. 4b), defined as [60]: For sufficiently long channels, the branch carrying most of discharge (say channel 2, continuous lines in Fig. 4) is subject to a gradient advantage with respect to the penalized branch ( S 2 > S 1 ); thus, discharge and slope asymmetry are concordant in sign. The model suggests that when the bifurcation parameter R is lower than the critical value R C a greater fraction of water and sediment fluxes is steered towards one of the two downstream branches as tidal amplitude increases (top half of Fig. 4a and b). In other words, tides act as a positive feedback sustaining the development of an unbalanced configuration. However, the scenario radically changes when the channels are shorter than a threshold length, L T . When L < L T , the signs of S 12 and Q 12 are discordant. In these conditions, tides exert a negative feedback by increasing the stability of the system, or rather, the possibility for the , respectively. The parameter R ∞ C denotes the critical threshold in the absence of tide. In panel b, the slope asymmetry index S 12 is plotted against the scaled channel length for a given value of the bifurcation parameter ( R = 0.2 ). The stabilizing/destabilizing effect of the tidal action for lengths shorter/longer than the threshold L T is highlighted for the case Q 12 < 0 . Yellow stars indicate threshold R C . A value of 0 = 20 is employed 1 3 system to partition fluxes evenly for lower values of R due to a steepening of the penalized branch (lower half of Figs. 4a and b).
It is worth observing that model results are similar if parameter 0 (see 3) instead of L is employed to analyze the response of the system. In general, for a given tidal amplitude > 0 , parameters 0 and L control how tidal action affects flow and sediment distribution through variations of the channel slopes, which are tracked by the asymmetry index S 12 . Such slope asymmetry is generated by the different values of the relative tidal amplitudes 1,2 at the internal nodes N 1,2 (see Eq. 11) that depend from how the fluxes are diverted upstream, in turn affecting the response of the bifurcation from downstream. In Fig. 5, the degree of asymmetry between the values of tidal amplitudes at the internal nodes 12 = ( 1 − 2 )∕ with respect to the configuration Q 12 < 0 (i.e. dominant channel 2) is shown. Interestingly, the most carrying branch is subject to a higher relative tidal amplitude (i.e. 2 > 1 ), regardless of the channels length being longer or shorter than the threshold length L T . Tidal asymmetry displays a minimum for a value of length L > L T , which roughly corresponds to the peak of the slope asymmetry index showed in Fig. 5. It follows that the tidal effect on the stability of the steady-state solutions is governed by the sign and magnitude of the slope asymmetry index S 12 . Figure 6 shows bed and water surface elevation profiles for channels 1 and 2 computed by means of Eqs. (8a-b), again for Q 12 < 0 . In the absence of tides (dashed lines), both channels have the same slope and water level but show a clear bed differential, with the dominating branch characterized by the larger depth. As expected, model outcomes coincide with the sole bifurcation problem as obtained by BRT [9]. Differently, bed and water surface elevations are a function of the position along the channels when > 0 , with results showing a substantial dependence on the scaled length being longer (top half of Fig. 6) or shorter (lower half of Fig. 6) than the threshold length L T . In the latter case, the penalized branch (i.e. channel 1) undergoes a progressive steepening moving downwards, which is induced by a higher tide-induced scour, till eventually becoming deeper than the most-carrying branch at the node (i.e. 1 | N 1 > 2 | N 2 ). On Let us now analyze what happens in the terminal branches (i.e. 3 to 6). In this case, the channels are constrained by the same value of the tidal amplitude at the mouths. As such, the effect of tide is always stabilizing as found in the single bifurcation case by RTB (Fig. 7) [57]. The value of the threshold at which the bifurcation evolves toward an unbalanced equilibrium state increases (Fig. 7a), suggesting the tendency of the system to attain stable balanced configurations for lower values of R and higher tidal amplitudes. Differently from the apex node, the channel slope is found to be always higher in the penalized branch (Fig. 7b). In accordance with previous studies [33,35,57,66], an increasing tidal influence results in a more even discharge partitioning in the terminal branches (i.e. reduced values of the discharge asymmetry Q).
A note of caution should be cast in the interpretation of model results presented in Fig. 7. Indeed, flux distribution at the apex node highly affects flow conditions of seaward channels. Specifically, as shown in Fig. 8a when the flux distribution is highly unbalanced the sediment transport capacity in the lowest carrying branches decreases till eventually vanishing (e.g., channel 4). As also suggested by BRT, under these conditions a progressive infilling with a channel abandonment is highly expectable. Moreover, low values of the conveyed water and sediment fluxes make channels more sensitive to tidal influence. From the definition as given in (3), we can define a "local" tidal amplitude ̃= ∕D u quantifying the tidal magnitude in the seaward channels that is effectively felt. Figure 8b shows that despite = 0.1 , the shallowest and lowest carrying branch is characterized by a value ̃4 ≃ 1 suggesting that in such conditions it might be substantially controlled by tidal motion.

Discussion
The present analysis shows that a contrasting response of channel bifurcations to tidal action arises even in an idealized steady model of a delta network. Specifically, results presented in Sect. 4 suggest that there is a threshold of the reference parameters (i.e. R, 0 , , L ) that separates two drastically different responses when multiple coupled bifurcations are affected by tidal motion. For a given value of parameters 0 and R requiring that > 0 , if L is longer than a threshold L T tide promotes the allocation of a larger fraction of the incoming discharge towards the most-carrying channel. The development of an unbalanced flow distribution is physically sustained through variations of the channel slopes, which in turn modify flow rating curves and sediment transport capacity of the two bifurcates. On the other hand, a stabilizing effect is found for L < L T . In this case, the lowest-carrying branch is steeper than the dominant channel due to a negative feedback exerted by tides that tends to re-establish a balanced flow and sediment distribution. A spontaneous question is what values do the relevant dimensionless parameters controlling the stability of tide-influenced bifurcations typically attain in real deltas. Estimated values of ( R, 0 , L ) for several natural and experimental deltas are reported in Fig. 9 as a function of the scaled tidal amplitude . Specifically, the analysis is based on two sets of laboratory experiments, performed by Lentsch et al. [42], and sixteen natural rivers forming deltaic distributary networks, which are characterized by different climatic settings (e.g., the Arctic Yukon, the Mediterranean Nile), morphology (e.g., lobe size and shape, number of distributary channels), and marine processes (i.e. waves and tides). We note that deltas are selected with the only purpose of understanding the typical range of magnitude of the different parameters. In this respect, we do not rely on the classic triad of river-, wave-, or tide-dominance [27,36,46]. Nevertheless, we avoid to consider deltas where the value of the scaled tidal amplitude is of O(1) or higher, as they exceed the limit of validity of the RTB model. Values of the parameter 0 spans within the range 10 0 − 10 2 (Fig. 9a), whereas R ranges between 10 −1 and 10 −3 (Fig. 9b). Finally, as a reference length scale of the maximum distance where the apex node is placed, we investigate the order of magnitude of the so-called "avulsion length" ( L A ∶= L * A ∕L * B0 ). The avulsion length, defined as the distance from the coastline of the apex node [16,36], varies between 10 −1 and 10 0 (Fig. 9c). In the light of the present analysis, the range of values attained by the reference parameters would suggest that even micro-tidal conditions are sufficient to substantially affect the overall delta dynamics in the long-term.
It is worth noting that as the majority of deltas analyzed in Fig. 9 are constituted by gentle-slope fine-grained rivers characterized by high-mobility conditions (i.e. high values of the Shields number), the low values of R are not surprising (see Eq. 1). According to the sole bifurcation BRT model, we should expect that almost all the bifurcations in the analyzed deltas are strongly unbalanced. Consequently, the uneven distribution of flow and sediment transport should reflect in a substantial asymmetry of the channel width in the distributaries of the networks. Nonetheless, this is not the case when tidal forcing is accounted for. Considering for example the case of the Mahakam Delta, a number of field and numerical studies have shown how bifurcations in the network are characterized 1 3 by a more uniform discharge distribution than the one that would occur in a purely riverine environment [33,68]. Tide distinctly shows its fingerprint on the delta morphology with funnel-shaped estuarine channels in the seaward reaches often disconnected from the main distributaries [68]. Despite the strongly idealized character of the present theoretical framework prevents a direct comparison with real-world systems, model predictions reproduce, at least qualitatively, some of the observed features. Some instructive indications can be gathered by analyzing the values of the reference parameters. Specifically, despite an extremely small value of the bifurcation parameter R (i.e. R ≃ 0.03 ), the magnitude of the parameters tied to tidal action, namely ≃ 0.2 , 0 ≃ 85 , and L A ≃ 0.1 , would suggest that tide is sufficiently strong to actively affect the upstream node. Computed values of the critical parameters from the RTB model give an avulsion length L A that is shorter than the threshold length L T . Physically, the model indicates that the observed small inequalities in discharge partitioning are downstream controlled by the tidal influence, where the stabilizing effect of the bedload steering at the bifurcation has a negligible importance.
Furthermore, in a series of recent laboratory experiments Lentsch et al. [42] analyzed the effects of tides and sea-level rise on deltaic morphology. In pure riverine conditions (i.e. = 0 ), a widespread number of channels are observed to dissect the subaerial surface of the delta where frequent avulsion events are triggered by channels aggradation; the shifting of the maximum sediment flux among branches leads over time to a nearly symmetric planform shape of the delta. Differently, as tide is progressively increased the number  [11,42,48,55,66,76,81] and lateral mobility of channels drop and a main thread is observed to distribute flow and sediment among few smaller branches originating towards the coastline. In this case, from the model analysis we obtain a value of the bifurcation parameter R below the critical threshold R C ; this suggests that in the absence of tides any bifurcation should evolve asymmetrically. Moreover, in spite of the moderate tidal amplitudes employed in the laboratory series ( 0.4 < < 1 ) the computed values of 0 ≃ 2 and L A ≃ 1 indicate that tidal influence on bifurcation dynamics is rather modest; in these conditions L > L T and the development of an unbalanced configuration at the apex node is potentially sustained by tidal action. Several field [33,49,76], experimental [42] and numerical [28,62] studies showed that tide-influenced deltas are often characterized by an initially-formed main (or few) active stem from which smaller distributary channels branch out, where tide acts to distribute uniformly flow and sediment across the network. The scouring action performed by tidal currents prevents bed aggradation, the main responsible for the formation/reoccupation of new/old channels through avulsion events in river-dominated deltas with a negligible tidal influence, which are often dissected by a widespread number of branches radiating from the apex bifurcation [61] and just few channels transport the largest amount of water and sediment [22].
While a systematic analysis is beyond the feasibility of the present work, the relative simplicity of the proposed theoretical model lends itself to be tested and validated through appropriate experimental and numerical analysis. In particular, those analysis would allow to obtain evidence of the two-way coupling between bifurcations in a deltaic network disentangling the role of allogenic mechanisms in the discharge-partitioning dynamics. As suggested by [65], an interaction between bifurcations could be also expected if the channels are allowed to establish their width freely or the delta geometry is allowed to change through time.
Let us discuss the main limitations of the present analysis. First, the potential effect of short-term temporal changes as daily fluctuations in the water levels, different tidal propagation in the branches [30,80,82], or due to wind-driven currents [70], are not accounted for by the model, which focuses on the long-term morphodynamic equilibrium state towards the system would tend. From a modeling perspective, relaxing the steady assumption would lead to some interesting developments, which are inspired from field and numerical evidences of the deltas discussed above. For example, during a sequence of spring-neap cycles in a couple of bifurcations in the Mahakam Delta the discharge division was observed to oscillate around an asymmetric state [66,67]. Differently, in the case of Wax Lake Delta, which represents a benchmark for the study of human-based solutions aimed at counteracting the sinking of many deltaic coastlines due to climate change, subsidence, and sediment starvation [14,15,32,51], during spring-neap cycles the discharge partitioning over two branches was observed to oscillate around a nearly symmetric state (see Fig. 4 of [80]). An unsteady analysis would allow to analyze the response of the system to time-varying forcing in a similar fashion to the work of Bertoldi et al. [7]. Bertoldi and co-workers studied experimentally and theoretically the effect of periodic downstream-migrating alternate bars interacting with a single fluvial bifurcation, which are shown to enrich and complicate the flux-partitioning dynamics. In a complementary way, throughout the paper we have recalled how flow oscillations can be generated even without time-dependent external forcing, but rather due to differential deposition in the branches [21,63]. A natural extension of the present model should include sediment deposition, incorporating into the single bifurcation model of Salter et al. [63] the tide-averaged effect following the framework of RTB. Second, we have employed a total load predictor for sediment transport. A natural extension would require to divide the material load into a bedload and suspended load component and then account for the gravitational effect on the lateral transport just on the former fraction.
Third, often tide-influenced rivers present a funnel-shaped morphology characterized by a gradual width expansion moving seaward [17,18,40]. Therefore, the present model may underestimate the potential destabilizing role of channel widening, which leads to increasing values of the width-to-depth ratio (i.e. lower values of R ) associated to channel shallowing. In a recent study, the drop of depth induced by the flow expansion has been observed to typically overwhelm the tide-induced scour [58]. Furthermore, we limited our analysis to the free response of the system, thus considering a symmetrical configuration where the connecting channels have the same width and length. However, in nature the geometry of the branches in a network is seldom the same. Different lengths affect the slope asymmetry and, as a consequence, the flow rating curves and sediment transport in the channels. Differently, in the present work variations of the channels slope arises internally in a purely symmetric configuration due to a mutual response in the form of a downstream-upstream feedback between the bifurcation and the confluence.
Fourth, the present formulation based on the two-cell model of BRT does not take into account the possible role of bifurcation angle. However, as shown by the two-dimensional linear analysis performed by [59], flow partitioning upstream the bifurcation seems to be poorly affected by an angle-induced effect. Indeed, the main physical mechanism generated by the angle is the formation of a separation zone induced by the sharp deviation of the flow at the node. However, it could be expected that this effect would be important for the local morphology in the downstream channels.
Finally, along the paper we have avoided to discuss what happens in tide-dominated conditions (i.e. ≫ O(1) ), representative of deltas like the Ganges-Brahmaputra-Meghna Delta, the Indus Delta, the Niger Delta [4,17,54], or in the coastal fringe of the Mahakam Delta [33]. Such networks diverge from the tree-like structure sketched in Fig. 3. The tidal plain is often characterized by tidal channels growing landward from the shoreline that frequently merge into secondary distributaries carrying freshwater and sediment, in turn generating looping patterns [39]. In these cases, the framework discussed herein must be abandoned in favor of a scheme explicitly including the interplay between riverine and bi-directional tidal flow. We note that a bifurcation model in purely tide-dominated conditions (i.e. in the absence of a fluvial supply of freshwater and sediment) is still lacking.

Conclusions
In this work, we have explored the long-term equilibrium configuration of a tide-influenced delta through a one-dimensional theoretical modeling. Though the theoretical model has been formulated for a tidal network characterized by an arbitrary bifurcation order k, in order to understand the basic mechanisms underlying the response of the system we have limited our analysis to a simple network composed by two orders of bifurcations. On the basis of model results, the following key outcomes are drawn: -Differently from the single bifurcation case, tides can be either a stabilizing or a destabilizing factor for the asymptotic equilibrium state reached by the system. The model indicates that this contrasting response crucially depends on slope variations in the channels induced by the different relative tidal amplitude in the internal nodes, which is a function of the flow and sediment delivered by the apex bifurcation and the ratio between the length of the channels and the backwater length; -The pluralistic behavior showed by the upstream node is driven by the sign of the slope asymmetry S with respect to the discharge asymmetry Q . In particular, if the scaled length of the channels is longer than a threshold length L T , S and Q are concordant in sign, and an initially stable symmetric bifurcation tends to become increasingly biased towards an unbalanced configuration as the amplitude of tidal oscillations is increased. This behavior is associated to an increased tidal strength in the branch receiving more discharge, which tends to steepen with respect to the other branch. Nonetheless, if the tidal influence is sufficiently strong and L < L T , S changes in sign. Consequently, the response of the system reverses, with tide hindering the abandonment of one of the two anabranches by increasing the slope of the lowest-carrying branch, i.e. providing a stabilizing effect; -Tides always exert a negative feedback in the terminal branches, reducing the inequalities in discharge partitioning that would occur in a pure riverine case; -In accordance with previous works [53,65], the model shows that tide-influenced delta networks behave as complex systems [45], where the overall response cannot be inferred without taking into account the coupling between upstream and downstream bifurcations.

Appendix A: Analytical expression of the coefficients , Ψ , and 2
In this section we provide the complete expressions for the coefficients appearing in Eqs. (8a-b). These coefficients are a function of the basic flow through the parameter 0 , and of the specific channel characteristics through the dependence on D u and S: where i is the imaginary unit, and the functional dependence of = f ( 0 , D u , S∕S 0 ) is given in the main text. We note that (15a-c) contribute to the magnitude of the tidal second-order correction (i.e. ∝ 2 ) that deviate flow variables (i.e. D, , U ) from the uniform conditions.

Appendix B: Tidal amplitude in the internal nodes
The starting point for the derivation of Eq. (11) is the one-dimensional differential problem for the evolution of a single river-dominated estuary. Namely, the governing equations for each branch of the network are the continuity and momentum conservation for the fluid phase, coupled with the sediment mass balance: where q and qs are the flow and sediment fluxes per unit width, p is the sediment porosity, x = x D −1 u S∕S 0 is a longitudinal coordinate with the origin at the sea (or at an internal node) and pointing landward, and T = q * 0 ∕ √ g * d * 3 is the ratio between the hydrodynamic timescale and "morphodynamic" timescale.
In the assumptions of micro-tidal conditions (i.e. ≪ 1 ) and F 2 0 ∼ F , with F an order O(1) quantity, Seminara et al. [69] tackled system (16a-c) to find an analytical solution for the tidally-averaged morphodynamic equilibrium by means of a perturbation approach. We recall that the equilibrium is achieved if: with the brackets denoting a tidal averaging. Essentially, Eq. (17) states that the net sediment flux in a tidal cycle is constant and coinciding with the sediment flux delivered by the channel.
Following the work of [69], each variable of the problem can be expanded in powers of as follows: where we note that now the suffixes denote the subsequent perturbation orders. Substituting from (18) into the governing system (16a-c), equating likewise powers of , a cascade of differential problems is then found. Here, with the aim of finding the superelevation of the free surface with respect to the reference level (i.e. H 0 ∝ L S∕S 0 ) moving landward along a channel from the shoreline, we are just interested in the order O( ) part of the solution, which describes the propagation of a small amplitude wave. For the interested readers, we refer to [69] for the details of the analysis. The solution is expressed as: where t = t * * , and c.c. stands for the complex conjugate of a complex number. Note that the tidally-averaged discharge per unit width q does not differ from the leading-order value q 0 , thus q 10 (i.e. the steady component) is zero.