2.1 The theory of morphodynamic influence

Two-dimensional morphodynamic influence can be defined as the process whereby the presence of a local persistent perturbation in a channel geometry is felt downstream and/or upstream through a series of damped, steady bed oscillations typically taking the form of alternate bars. Theoretically, this emerges from the linear theory of non-migrating bars in straight channels with a non-uniform boundary condition at one end of the channel, which translates in mathematical terms to a discontinuity. This phenomenon has been known for decades to occur downstream of the local perturbation (Struiksma et al. Reference Struiksma, Olesen, Flokstra and de Vriend1985; Struiksma & Crosato Reference Struiksma, Crosato, Ikeda and Parker1989), while the discovery that, under suitable conditions, morphodynamic influence can also be felt upstream of the geometrical perturbation is more recent (Zolezzi & Seminara Reference Zolezzi and Seminara2001). Such a phenomenon has been experimentally verified (Zolezzi et al. Reference Zolezzi, Guala, Termini and Seminara2005) and numerically reproduced (van der Meer et al. Reference van der Meer, Mosselman, Sloff, Jager, Zolezzi and Tubino2011; Siviglia et al. Reference Siviglia, Stecca, Vanzo, Zolezzi, Toro and Tubino2013), and has been mostly associated with a localised persistent perturbation, such as a local narrowing in a straight channel or a discontinuity in channel curvature in a stream of constant width.

A crucial parameter of the theory is the channel aspect ratio, ${\it\beta}$ , defined as the half-width-to-depth ratio. It is now broadly agreed (Mosselman et al. Reference Mosselman, Tubino and Zolezzi2006) that 2-D morphodynamic influence can occur upstream of the localised perturbation only if ${\it\beta}$ is larger than a threshold resonant value ${\it\beta}_{R}$ , whereas only downstream influence occurs when ${\it\beta}<{\it\beta}_{R}$ . The term ‘resonant’ refers to the unbounded response of the linear bar solution in meandering channels with periodic distribution of curvature as ${\it\beta}$ matches ${\it\beta}_{R}$ (Blondeaux & Seminara Reference Blondeaux and Seminara1985).

The legacy between morphodynamic influence and bifurcation dynamics can be understood when viewing the bifurcation node as a localised persistent perturbation found at the downstream end of the bifurcating channel and at the upstream ends of the two downstream bifurcates.

2.2 Channel bifurcations: simplified modelling and morphodynamic influence

Modelling of channel bifurcation morphodynamics has been proposed in the last three decades through simplified one-dimensional (1-D) models for each of the three distributaries, plus suitable boundary conditions to be imposed at the upstream end of the bifurcating channel, at the downstream ends of the bifurcates and at the bifurcation node. Namely, adopting a 1-D mobile-bed model for every channel requires us to specify two conditions at the upstream end of each distributary and one condition at the downstream end of the inlet channel, for a total of five nodal relations (Bolla Pittaluga et al. Reference Bolla Pittaluga, Repetto and Tubino2003). In addition to the conservation of water and sediment mass, most 1-D models (e.g. Bolla Pittaluga et al. Reference Bolla Pittaluga, Repetto and Tubino2003; Miori, Repetto & Tubino Reference Miori, Repetto and Tubino2006; Kleinhans et al. Reference Kleinhans, Cohen, Hoekstra and Ijmker2011; Bolla Pittaluga et al. Reference Bolla Pittaluga, Coco and Kleinhans2015) impose the continuity of free surface elevation at the node. Setting the fifth relation is more complex because it needs to take into account the structure of the flow field and sediment transport near the bifurcation node. The type of this boundary condition and the values of the involved parameter are key to discriminating between contrasting types of bifurcation behaviours. Moreover, the nodal condition is nearly the only opportunity of incorporating 2-D information related to water and sediment partition at the bifurcation within a 1-D approach.

The nodal condition proposed by Bolla Pittaluga et al. (Reference Bolla Pittaluga, Repetto and Tubino2003) provides a simple description of the 2-D processes at the bifurcation node based on the definition of two cells (figure 1) of length ( ${\it\alpha}W_{A}^{\ast }$ ) proportional to the width of the inlet channel, which exchange water and sediment fluxes depending on their bottom elevation differences.

Figure 1. Sketch of a channel bifurcation. Cells 1 and 2 refer to the Bolla Pittaluga et al. (Reference Bolla Pittaluga, Repetto and Tubino2003) model.

The model indicates that, when the aspect ratio of the inlet channel ${\it\beta}_{A}$ is higher than a critical threshold, whose value depends on the flow and sediment parameters of the incoming flux, the balanced solution is unstable and an uneven configuration arises. Such an unbalanced configuration has two key characteristics: the discharge partition between the two bifurcates is uneven and there is an elevation gap (‘inlet step’) between the bed levels at the inlet of the two bifurcates.

This simple model is able to capture important bifurcation properties observed both in laboratory physical models and in the field (Ferguson et al. Reference Ferguson, Ashmore, Ashworth, Paola and Prestegaard1992; Zolezzi et al. Reference Zolezzi, Bertoldi, Tubino, Best, Bristow and Petts2006). It predicts the development of unbalanced bifurcation configurations for relatively low values of Shields stress and large aspect ratio, as found in many individual branches of gravel-bed bifurcations.

The main disadvantage of this quasi-2-D approach is that the resulting critical aspect ratio is highly dependent (i.e. directly proportional) on the parameter ${\it\alpha}$ , which measures the upstream reach length affected by the bifurcation. Calibration of ${\it\alpha}$ often yields rather different values: $1$ in Bolla Pittaluga et al. (Reference Bolla Pittaluga, Repetto and Tubino2003), $2.5$ in Siviglia et al. (Reference Siviglia, Stecca, Vanzo, Zolezzi, Toro and Tubino2013), $3$ in Kleinhans et al. (Reference Kleinhans, Jagers, Mosselman and Sloff2008) and $6$ in Bertoldi & Tubino (Reference Bertoldi and Tubino2007).

The clearest and perhaps most intriguing indication of the factors controlling the length of the bifurcation-affected upstream channel reach follows from the experimental analysis of Bertoldi & Tubino (Reference Bertoldi and Tubino2007). They observed that the value of the parameter ${\it\alpha}$ which yields the best matching with the measured discharge partition is highly correlated with the difference between the bar-forming value of the inlet channel aspect ratio ${\it\beta}_{A}$ and the resonant threshold ${\it\beta}_{R}$ . Such observation sets the connection between the theory of morphodynamic influence and the dynamics of bifurcations. In particular, as reported in figure $12$ of Bertoldi & Tubino (Reference Bertoldi and Tubino2007), the optimal cell size ( ${\it\alpha}W_{A}^{\ast }$ ) shows its peak under near-resonant conditions; this is consistent with the theory which predicts the spatial damping of the steady bars to vanish at the resonant point. In addition Bertoldi & Tubino (Reference Bertoldi and Tubino2007) noticed that the final configuration reached at the end of experimental runs crucially depends on the distance from the resonant conditions; indeed, the two indicators of bifurcation asymmetry, namely the discharge ratio and the difference in bed level at the inlet of downstream anabranches, show a strong correlation with the scaled difference $({\it\beta}_{A}-{\it\beta}_{R})/{\it\beta}_{R}$ rather than with the aspect ratio itself.

3.1 The steady two-dimensional model

We adopt a Cartesian system of coordinates $\{x^{\ast },y^{\ast }\}$ and we call $U^{\ast }$ and $V^{\ast }$ the longitudinal and cross components of the depth-averaged velocity vector (figure 2 a). We define $\{{\it\tau}_{x}^{\ast },{\it\tau}_{y}^{\ast }\}$ and $\{\boldsymbol{q}s_{x}^{\ast },\boldsymbol{q}s_{y}^{\ast }\}$ the components of the bed shear stress and of the specific (per unit width) solid discharge, respectively. In addition we use the symbols $D^{\ast }$ for water depth and ${\it\eta}^{\ast }$ for bed elevation (figure 2 b).

The above quantities can be made dimensionless as follows:

where ${\it\rho}$ is the water density, $B^{\ast }$ is the horizontal length scale, taken to coincide with the width of the bifurcates, $U_{0}^{\ast }$ , $D_{0}^{\ast }$ and $qs_{0}^{\ast }$ are reference values for velocity, depth and unit sediment discharge corresponding to uniform flow conditions, for given discharge, sediment size and average slope.

In the following we investigate the steady bifurcation configuration: therefore we neglect the time derivatives and we focus on the spatial structure of the solution. The 2-D, depth-averaged conservation of $x$ -momentum, $y$ -momentum, water mass and sediment mass can be expressed through the following set of partial differential equations:

where the Froude number and the aspect ratio are given by:

In order to close the problem we need to specify the bed shear stress and sediment transport components in terms of the primitive variables ${\it\eta},U,V,D$ . Following a widely adopted procedure, this is achieved through the ‘local equilibrium’ approximation, whereby closure relationships for the near-bed shear stress and bedload rates are evaluated in terms of the local values of the primitive variables. Specifically, the bed shear stress is expressed by the Chèzy formula and is assumed to be oriented as the velocity vector, namely

where the dimensionless Chèzy coefficient can be estimated as in Engelund & Fredsoe (Reference Engelund and Fredsoe1982):

and $d_{s}^{\ast }$ is the grain size.

Assuming bedload as the dominant mode of sediment transport we use the Einstein (Reference Einstein1950) scaling, namely

where ${\it\Delta}$ is the relative submerged density of the bed material and the Shields parameter ${\it\theta}$ is defined as

The dimensionless sediment flux is consequently expressed as

where ${\it\Phi}({\it\theta})$ is computed using a suitable bedload formula such as Parker (Reference Parker1990) and ${\it\theta}_{0}$ is the Shields stress of the reference uniform flow.

The direction of the sediment transport vector is given by:

with

where ${\it\gamma}_{q}$ represents the direction of the velocity vector, ${\it\gamma}_{g}$ is a correction due to the transverse bed slope, $n$ is the coordinate orthogonal to the local velocity vector and $r$ is an empirical constant typically ranging from $0.3$ to $0.6$ (Ikeda, Parker & Sawai Reference Ikeda, Parker and Sawai1982; Talmon, Van Mierlo & Struiskma Reference Talmon, Van Mierlo and Struiskma1995).

We note that a unique definition of the reference flow is given in dimensionless form once the independent parameters ${\it\beta}$ , ${\it\theta}_{0}$ , $d_{s}$ and ${\it\Delta}$ are fixed, as the following conditions hold for uniform flow conditions:

where

and $S_{0}$ is the average longitudinal slope.

3.2 The steady linear solution for a straight channel

The governing dimensionless mathematical model consists of a system of four nonlinear partial differential equations in two space dimensions (3.2) whose initial boundary problem can be fully solved only through numerical approximation (e.g. Siviglia et al. Reference Siviglia, Stecca, Vanzo, Zolezzi, Toro and Tubino2013). Nevertheless, under the hypothesis of relatively small perturbations with respect to a known basic solution, it is possible to linearize the problem, which enables us to obtain an analytical solution for some particular geometrical and boundary configurations. The analytical approach makes it easy to detect the key controlling parameters and to efficiently explore the system behaviour within a broad range of controlling factors.

The basic solution is taken as the uniform flow that would occur on a flat-bottom sloping channel, namely:

A perturbation solution of (3.2) is obtained by expanding the four unknowns $({\it\eta},U,V,D)$ around the basic state as follows:

where the small parameter ${\it\epsilon}$ measures the order of magnitude of the perturbations and $O({\it\epsilon}^{2})$ indicates quadratic and higher-order terms that can be neglected at a first approximation as far as perturbations are sufficiently small to ensure that nonlinear effects are weak.

The complete linear solution for the free bifurcation problem requires us to separately solve the 2-D linear morphodynamic model for each of the three channels $A$ , $B$ , $C$ and then to match the solutions through suitable nodal conditions. Matching also implies that every solution is obtained referring to the same Cartesian coordinate system ( $x,y$ ) displayed in figure 2(a). This means that the origin of the $y$ -axis is located on the centreline of channel $A$ , on the right bank of channel $B$ and on the left bank of channel $C$ .

In this section we review the steady linear solution as obtained by Zolezzi & Seminara (Reference Zolezzi and Seminara2001) with reference to the main channel  $A$ . Once the solution for channel $A$ is obtained, calculating the solution for the sub-domains $B$ and $C$ is straightforward because it only requires us to shift the location of the streamwise axis $x$ with respect to the lateral channel boundaries (see §§ 4.1 and 4.2 for more details).

The linear approximation of the system (3.2), at the order $O({\it\epsilon})$ , reads:

where ${\it\Phi}_{T}$ and $C_{D}$ express the dependence of bedload and flow resistance on Shields stress and depth, respectively, and read:

The first-order linear system (3.15) can be solved by Fourier analysing the unknowns in the transverse ( $y$ ) direction, thus obtaining a cascade of fourth-order ordinary differential problems in the independent variable ( $x$ ) for each Fourier mode $m$ . Linearity of the system enables us to solve the differential problem for each mode $m$ separately.

Following Zolezzi & Seminara (Reference Zolezzi and Seminara2001) we then write the general solution for each Fourier mode $m$ in the form:

where $\text{c.c.}$ denotes the complex conjugate and the longitudinal variation is defined as

In (3.18) ${\it\lambda}_{mj}$ is a complex wavenumber, whose real part represents the spatial growth rate of perturbations, while the imaginary part defines their spatial periodicity. Furthermore, the transverse structure of the different modes reads:

We note that, unlike Zolezzi & Seminara (Reference Zolezzi and Seminara2001), who considered only the odd modes (the asymmetrical component of the solution), we need in general to include also the even modes (symmetrical component), as required to represent the bifurcation configuration.

Substituting (3.17)–(3.19) into the linear system (3.15) we obtain a fourth-order polynomial of the form

which in general admits four complex eigenvalues ${\it\lambda}_{mj}$ $\left(j=1,\ldots ,4\right)$ , for each Fourier mode $m$ .

The coefficients $\{u_{mj},v_{mj},d_{mj}\}$ can be computed from system (3.15) in terms of the parameters of the basic flow, while the complex constants $\tilde{{\it\eta}}_{mj}$ , being undetermined, can be arbitrary chosen.

In the following we will derive a solution in the complex field and we will therefore omit the complex conjugate notation for the sake of clarity. However, we should keep in mind that only the real part of the final solution for the variables $\{{\it\eta}_{1},U_{1},V_{1},D_{1}\}$ is physically significant.

3.3 Solution of the linear system and morphodynamic influence

As pointed out by Zolezzi & Seminara (Reference Zolezzi and Seminara2001), the sign of the real part of ${\it\lambda}_{j}$ (spatial growth rate) of each solution plays a key role in determining the dominant direction of the morphodynamic influence that can be induced by the presence of a local discontinuity. In fact, spatially growing perturbations that do not vanish far from the discontinuity are not compatible with the conditions set by the reference flow in a semi-infinite domain.

The dispersion relation (3.20) implies that the values of complex eigenvalues change with the whole set of parameters of the reference flow, with the spatial growth rate being mainly controlled by the aspect ratio ${\it\beta}$ , for given values of ${\it\theta}_{0}$ and $d_{s}$ . The example reported in figure 3 illustrates the dependency on the channel aspect ratio ${\it\beta}$ of the four eigenvalues ( ${\it\lambda}_{1}$ , ${\it\lambda}_{2}$ , ${\it\lambda}_{3}$ , ${\it\lambda}_{4}$ ) of the 2-D steady solution: results are given with reference to the first transverse mode only ( $m=1$ ), therefore the index $m$ is omitted to simplify the notation. Within a wide range of ${\it\beta}$ (say 5–18), two eigenvalues ( ${\it\lambda}_{1}$ and ${\it\lambda}_{4}$ ) are purely real and relatively large, while the other two eigenvalues ( ${\it\lambda}_{2}$ and ${\it\lambda}_{3}$ ) are complex conjugates and exhibit a comparatively smaller growth rate. The two real eigenvalues represent a local effect which grows (or decays) rapidly in $x$ direction, while the two complex eigenvalues embody the fundamental 2-D morphodynamic influence which manifests itself through the formation of a steady, spatially growing (or damped), alternate bar pattern. The above scenario does not change qualitatively when different values of ${\it\theta}_{0}$ and $d_{s}$ are used.

Figure 3. Dependence on the channel aspect ratio ${\it\beta}$ of the four eigenvalues ${\it\lambda}_{mj}$ of the eigenrelation (3.20) for the first transverse mode ( $m=1$ ) and $d_{s}=0.05$ , ${\it\theta}_{0}=0.1$ , $r=0.5$ : (a) real part (spatial growth rate); (b) imaginary part (wavenumber).

Figure 3(a) shows that the real part of the two complex conjugate eigenvalues $({\it\lambda}_{2})_{R}=({\it\lambda}_{3})_{R}$ changes sign from negative to positive as the aspect ratio ${\it\beta}$ crosses the resonant threshold ${\it\beta}_{R}$ . This condition sets the fundamental connection between resonance and morphodynamic influence, which turns out to be the key ingredient to interpreting the dynamics of the bifurcation. Two distinct ‘morphodynamic regimes’ can be identified:

1. (i) super-resonant ( ${\it\beta}>{\it\beta}_{R}$ ), dominant upstream influence, when the real part of three ( ${\it\lambda}_{2}$ , ${\it\lambda}_{3}$ and ${\it\lambda}_{4}$ ) out of four eigenvalues is positive;

2. (ii) sub-resonant ( ${\it\beta}<{\it\beta}_{R}$ ), dominant downstream influence, when only one eigenvalue ( ${\it\lambda}_{4}$ ) has positive real part.

To better clarify such a connection we first consider a semi-infinite rectangular channel extending to $-\infty$ (left panel in figure 4), which mimics channel $A$ joining at the node at the downstream end. If we look for bounded (i.e. limited) solutions set by boundary conditions in $x=0$ and having the form of (3.17), we must discard those solutions having $({\it\lambda}_{j})_{R}<0$ as they grow exponentially in the upstream (negative  $x$ ) direction. Consequently, in the sub-resonant regime ( ${\it\beta}<{\it\beta}_{R}$ ), the only acceptable solution, corresponding to the eigenvalue ${\it\lambda}_{4}$ , is rapidly damped in the upstream direction, and therefore the upstream influence is felt only locally. On the other hand, in the super-resonant regime ( ${\it\beta}>{\it\beta}_{R}$ ) two further independent solutions are acceptable, corresponding to the complex conjugate eigenvalues ${\it\lambda}_{2}$ and ${\it\lambda}_{3}$ , which implies that the influence of boundary conditions can be felt relatively far upstream through the formation of damped alternate bars.

The opposite behaviour characterizes a semi-infinite channel which extends towards $+\infty$ (right panel of figure 4), as in the case of channels $B$ and $C$ joining at the node at the upstream end. In this case the occurrence of a spatially damped alternate bar pattern in the downstream direction is only possible in the sub-resonant regime, where $({\it\lambda}_{2})_{R}=({\it\lambda}_{3})_{R}$ is negative, while in the super-resonant regime ( ${\it\beta}>{\it\beta}_{R}$ ) the only acceptable solution (associated with the eigenvalue ${\it\lambda}_{1}$ ) produces a downstream influence which vanishes quite rapidly.

Figure 4. Structure of the solution for semi-infinite channels under sub-resonant (a) and super-resonant (b) conditions. BC denotes the section ( $x=0$ ) at which the boundary conditions are applied, $m$ indicates the transverse Fourier mode and $j$ the eigenvalue.

The results reported in figure 3 can be readily extended to the solutions for higher ( $m>1$ ) modes when considering that the dispersion relation (3.20) obeys the following scaling property:

which indicates that the solution for the higher ( $m>1$ ) modes is identical to that for the first mode ( $m=1$ ) if we take ${\it\lambda}_{mj}=m{\it\lambda}_{1j}$ and scale the aspect ratio ${\it\beta}$ by a factor  $m$ .

Therefore, according to (3.21), in both sub-resonant and super-resonant regime, provided ${\it\beta}<2{\it\beta}_{R}$ , three eigenvalues have a negative real part for all $m>1$ , which means that for the higher transverse modes there are always three linearly independent solutions in the downstream semi-infinite channel and one in the upstream semi-infinite channel of figure 4. In other words, as long as ${\it\beta}$ does not exceed $2{\it\beta}_{R}$ , higher-order modes behave as in the sub-resonant regime and therefore the occurrence of steady $m>1$ bed oscillations is only possible in the downstream direction.

3.4 The one-dimensional component

The linear solution of system (3.15) also includes a 1-D component $\{{\it\eta}_{1},U_{1},D_{1}\}^{1D}$ corresponding to the transverse mode $m=0$ , which cannot be written in the form of (3.17).

Discarding the $y$ -derivatives in the continuity equations for water and sediment (3.15c ) and (3.15d ) one readily finds that the 1-D solution is uniform, namely

Substituting into the longitudinal momentum equation (3.15a ) and recalling (3.11b ) gives

where the (arbitrary) integration constant $\tilde{{\it\eta}}_{0}$ indicates a uniform variation of the bottom elevation and $S_{1}$ represents a perturbation of the bed slope. Therefore, the 1-D component is determined once the three constants $U_{1},D_{1},\tilde{{\it\eta}}_{0}$ are fixed.

In order to better highlight the physical meaning of such a component, it is convenient to re-write (3.22) and (3.23) as follows:

where the coefficients ${\it\gamma}_{1}$ and ${\it\gamma}_{2}$ account for the variations of velocity and depth due to a perturbation of the discharge ( $Q_{1}=U_{1}+D_{1}$ ), with fixed-bed elevation, and read

while the coefficient ${\it\gamma}_{3}$ measures the perturbation of depth and velocity due to slope variation, with fixed discharge, and is defined as

The 1-D component (3.24) is constant in the transverse direction and physically expresses the lowest-order ( $m=0$ ) channel response which can be driven by a uniform variation of bottom elevation $\tilde{{\it\eta}}_{0}$ or flow discharge $Q_{1}$ , as well as by a longitudinal slope adjustment (term proportional to $S_{1}$ ). We note that such a component, which was not included in previous theories for steady bars in straight channels (e.g. Struiksma & Crosato Reference Struiksma, Crosato, Ikeda and Parker1989; Seminara & Tubino Reference Seminara and Tubino1992), is relevant in the present case of a bifurcating configuration as it brings the fundamental ingredient in reproducing differences in discharge and bed level between the two distributaries.

The complete steady linear solution in a straight channel is finally obtained by summing the 2-D (3.17) and 1-D (3.24) components. Specifically, considering a finite sum of the first $N$ Fourier modes, the solution reads:

Equation (3.27) shows that the number of degrees of freedom for the steady linear solution in each channel, i.e. the number of complex constants that can be chosen arbitrarily, is equal to $4N-1$ , namely

4.1 Solution for the main channel A

The solution for the inlet channel  $A$ can be readily derived from (3.27) discarding all the even modes $m$ as in Zolezzi & Seminara (Reference Zolezzi and Seminara2001) and including the 1-D component. As illustrated in § 3.3, the structure of the solution depends on ${\it\beta}$ being smaller or greater than the resonant threshold ${\it\beta}_{R}$ . The two cases are then tackled separately in the following sub-sections.

4.2 Solution for the downstream branches B, C

The anti-symmetric structure of the problem implies that finding the solution in the bifurcates only requires us to consider one distributary (hereinafter the left channel  $B$ ), the corresponding solution in the other channel being the opposite.

Differently from the main channel  $A$ , the secondary channel  $B$ has half-width ( $B^{\ast }$ ) and the reference system is placed on the right bank (see figure 2 a). Therefore, the solution can be given the form (3.27), where

and the Fourier expansion includes both even and odd modes, as the solution in a single bifurcate does not follow any symmetry property.

Moreover, as $W_{A}^{\ast }=2W_{B}^{\ast }$ , the secondary branch is always under sub-resonant conditions, regardless the upstream channel flow being sub- or super-resonant, and provided ${\it\beta}_{A}$ does not exceed $2{\it\beta}_{R}$ . Therefore, as discussed in § 3.3, the exponentially growing solution corresponding to the fourth positive eigenvalue is not compatible with the semi-infinite character of the channel. Similarly, any change of slope $S_{1}$ produces a variation of the bottom level which becomes infinite when $x\rightarrow \infty$ , and is therefore not acceptable. Hence, the solution for the downstream bifurcate $B$ reads:

where the subscript $n$ is used to denote the transverse modes in downstream branches to avoid confusion with the $A$ -solution.

Solution (4.4) is a combination of $3N-1$ linearly independent components, whose amplitude is defined by the following constants:

4.3 The matching procedure

The general $N$ -modes solution derived for each channel is a combination of linearly independent components whose amplitude can be freely chosen. We now seek solutions which satisfy not only the differential problem within each channel, but also the matching condition across the boundary between inlet channel and distributaries.

As highlighted before, the general solution in the secondary channel  $B$ is a linear function of $3N-1$ independent parameters, while in the main channel  $A$ it is a combination of $N$ or $N+2$ , depending upon ${\it\beta}_{A}$ being smaller or larger than ${\it\beta}_{R}$ (see table 1). This difference in the number of degrees of freedom is key to understanding the different bifurcation behaviour under sub- and super-resonant conditions.

Matching is achieved by setting the $N$ -modes expansion of the main channel solution at $x=0$ equal to the solution in channel  $B$ , within the transverse range $y\in [0,1]$ . Due to the different width of the channels, the two solutions exhibit a different Fourier representation, as is apparent when comparing (4.1) or (4.2) with (4.4). This implies that the matching condition must be ensured by expanding each transverse $m$ mode of the upstream channel solution in terms of the Fourier representation used in the downstream $B$ channel (or vice versa).

For example, matching the longitudinal velocity at the node is imposed by expanding the upstream channel solution in the form:

where according to (4.1a ) and (4.2a ):

and the expansion coefficients are given by:

Equating (4.6) with the solution for the downstream bifurcates (4.4a ) gives:

where we have defined

Figure 5. Different Fourier representation of the solution between the main channel and the distributaries: each mode of the solution (4.2a ) for the longitudinal velocity in channel  $A$ (solid lines) is expanded using an increasing number $N$ of transverse modes of the solution in the downstream bifurcates (4.4a ) (dashed lines): (a) channel  $A$ mode $m=1$ ; (b) channel  $A$ mode $m=3$ .

The different Fourier representations of the solution between the main channel and the distributaries is highlighted in figure 5, where convergence to the exact solution for an increasing number of Fourier components is shown with reference to the first $m=1$ (a) and second $m=3$ (b) modes of the sine components of the $A$ -solution.

The matching procedure provides $N$ conditions of the form (4.9) for each of the variables ${\it\eta}$ , $U$ and $D$ , and $(N-1)$ conditions for the variable $V$ , as its solution in the downstream branches (4.4b ) does not contain any constant component $(n=0)$ .

Therefore, matching at $x=0$ results in an homogeneous system of $4N-1$ linearly independent equations. As the unknowns of this system are the ‘degrees of freedom’ reported in table 1, for ${\it\beta}<{\it\beta}_{R}$ , the number of equations equals the number of unknowns; consequently, the only solution is trivial and gives a vanishing perturbation. In other words, under sub-resonant conditions, the only steady solution compatible with the semi-infinite character of the channels is uniform flow. As illustrated in figure 6(a) the bifurcation node plays a ‘passive’ role and is therefore unable to produce a morphodynamic influence that propagates far upstream.

On the other hand, in the super-resonant case two more degrees of freedom (i.e. unknowns) exist, resulting in $4N+1$ unknowns, with $4N-1$ conditions. Therefore the system admits an infinite number of non-vanishing solutions, which are compatible with the conditions at $x\rightarrow \pm \infty$ . Formally, the homogeneous linear system exhibits a kernel of dimension two, so that $\infty ^{2}$ non-trivial solutions are possible; this implies that an unique solution can be found once two complex constants are independently fixed. Therefore, for ${\it\beta}_{A}>{\it\beta}_{R}$ , because upstream and downstream influences occur simultaneously (figure 6 b), the bifurcation node acts like a control section and can lead to non-vanishing unbalanced solutions.

4.4 The autogenic free response of bifurcations

The above results highlight the connection between bifurcation dynamics and morphodynamic influence. The bed elevation maps reported in figure 6 clearly show two distinct bifurcation configurations: (a) when ${\it\beta}_{A}<{\it\beta}_{R}$ (sub-resonant range) only the trivial uniform solution is possible, the bottom is flat and the discharge equally distributed; (b) when ${\it\beta}_{A}>{\it\beta}_{R}$ (super-resonant range) a pattern of steady damped alternate bars forms upstream of the bifurcation, which induces asymmetrical discharge partition and elevation difference (inlet step) between distributaries. Therefore, in the absence of external forcing, when ${\it\beta}_{A}<{\it\beta}_{R}$ a perturbation of flow and sediment partition at the node is not compatible with the steady linear solution and eventually vanishes, while for ${\it\beta}_{A}>{\it\beta}_{R}$ an unbalanced configuration arises as a result of a free, ‘autogenic’ response of the system. This is the main outcome of the present theory as it provides a rigorous proof to the experimental findings of Bertoldi & Tubino (Reference Bertoldi and Tubino2007), showing the intrinsic tendency of bifurcations to develop an asymmetrical configuration when the aspect ratio is large enough to enable upstream morphodynamic influence.

Figure 6. Direction of morphodynamic influence in the three channels joining at the node (arrows) and associated steady linear solution for bottom elevation ${\it\eta}$ (contour maps) for a bifurcation with vanishing angle (figure 2 a) and for two different conditions of the main channel $A$ : (a) sub-resonant ( ${\it\beta}_{A}=12$ ) case; (b) super-resonant ( ${\it\beta}_{A}=15$ ) case. The basic flow parameters are $d_{s}=0.02$ , ${\it\theta}_{0}=0.1$ , $r=0.5$ , the resonant threshold is ${\it\beta}_{R}=13.4$ .

Based on this result we can give a sound interpretation of the threshold value of the aspect ratio discriminating between balanced and unbalanced bifurcations as found by Bolla Pittaluga et al. (Reference Bolla Pittaluga, Repetto and Tubino2003), by referring to the resonant value ${\it\beta}_{R}$ . Plots of ${\it\beta}_{R}$ as a function of the Shields stress, ${\it\theta}_{0}$ , for different values of the dimensionless roughness $d_{s}$ , are reported in figure 4 of Seminara & Tubino (Reference Seminara and Tubino1992). In figure 7 we test the above criterion against the experimental data of Bertoldi & Tubino (Reference Bertoldi and Tubino2007), who classified ‘balanced’ or ‘unbalanced’ runs depending upon the measured discharge anomaly $(Q_{B}-Q_{C})/Q_{A}$ being smaller or larger than 0.03. Figure 7 also includes the results of 19 numerical runs performed by Siviglia et al. (Reference Siviglia, Stecca, Vanzo, Zolezzi, Toro and Tubino2013) with a fully nonlinear 2-D morphodynamical model, where the ‘balanced’ or ‘unbalanced’ character of the bifurcation was assessed by monitoring the temporal evolution of the discharge partition between the distributaries. In both cases the theoretical criterion performs reasonably well, as nearly the whole set of balanced data lie in the sub-resonant region, while the unbalanced data fall in the super-resonant domain, independently of the Shields stress.

Figure 7. Comparison between the theoretical criterion for balanced/unbalanced bifurcations based on the resonant aspect ratio ${\it\beta}_{R}$ and flume (Bertoldi & Tubino Reference Bertoldi and Tubino2007) and numerical (Siviglia et al. Reference Siviglia, Stecca, Vanzo, Zolezzi, Toro and Tubino2013) data.

Present theoretical findings are consistent with the observation that bifurcations in braided reaches of gravel-bed rivers, where super-resonant conditions are frequently met, are rarely symmetrical (Bertoldi & Tubino Reference Bertoldi and Tubino2007; Bertoldi Reference Bertoldi2012). Evidence of the relationship between morphodynamic influence and bifurcation dynamics can be obtained from field data reported by Zolezzi et al. (Reference Zolezzi, Bertoldi, Tubino, Best, Bristow and Petts2006) who observed the inherent (autogenic) tendency of the poorly vegetated, low-cohesion gravel-bed braided reaches of the Ridanna and Sunwapta rivers to produce unbalanced bifurcations, with uneven water distribution between distributaries, especially at relatively low flow stages. In table 2 the bankfull values of the aspect ratio of the incoming channel ${\it\beta}_{A}$ of the six unbalanced bifurcations analysed by Zolezzi et al. (Reference Zolezzi, Bertoldi, Tubino, Best, Bristow and Petts2006) are reported along with the corresponding computed resonant values ${\it\beta}_{R}$ . Sensitivity of the latter resonant values with respect to the chosen value of the empirical parameter $r$ , which accounts for gravitational effect on sediment transport (see 3.10b ), is also documented. We note that all bifurcations fall within the super-resonant regime, which suggests that their observed asymmetrical behaviour can be given a sound interpretation as a result of the capability of gravel-bed channels to undergo an upstream influence.

Figure 8. Values of the dimensionless length of upstream influence, ${\it\alpha}$ , of the Bolla Pittaluga et al. (Reference Bolla Pittaluga, Repetto and Tubino2003) model which best predict the discharge ratio observed in the Bertoldi & Tubino (Reference Bertoldi and Tubino2007) experiments (circle and square markers) are compared with the theoretical predictions, ${\it\alpha}^{T}$ of (4.11), based on the present theory (star markers).

As highlighted in § 2.2, further indication of the connection between bifurcation dynamics and the theory of morphodynamic influence was also given by Bertoldi & Tubino (Reference Bertoldi and Tubino2007) who computed the values of the coefficient ${\it\alpha}$ of the Bolla Pittaluga et al. (Reference Bolla Pittaluga, Repetto and Tubino2003) model which best predict the observed discharge ratios in the distributaries in their experiments. Being ${\it\alpha}$ in the theory of Bolla Pittaluga et al. (Reference Bolla Pittaluga, Repetto and Tubino2003) a measure of the upstream distance (scaled with $2B^{\ast }$ ) of the 2-D morphodynamical effect of the bifurcation, it should behave as the inverse of the damping rate $({\it\lambda}_{2})_{R}=({\it\lambda}_{3})_{R}$ , which vanishes at ${\it\beta}_{A}={\it\beta}_{R}$ . Therefore, ${\it\alpha}$ is expected to be infinite close to resonant conditions and to decrease as ( ${\it\beta}_{A}-{\it\beta}_{R}$ ) increases. In figure 8 computed values of ${\it\alpha}$ are compared with a theoretical estimate of such distance, ${\it\alpha}^{T}$ , set as the dimensionless length at which the amplitude of the linear solution (3.17) is reduced by one-third, namely

Note that ${\it\alpha}^{T}$ in figure 8 does not decrease monotonically due to the different values of $d_{s}$ and ${\it\theta}_{0}$ in the experimental runs.

4.5 The divergent symmetrical bifurcation

In the previous § 4.4 we have derived a linear analytical solution for the straight bifurcation geometry of figure 2, showing that non-vanishing perturbed solutions are only possible when ${\it\beta}_{A}>{\it\beta}_{R}$ . In this section we show that this scenario is not affected by the presence of a symmetric angle between the distributaries, as in the experimental runs of Bertoldi & Tubino (Reference Bertoldi and Tubino2007) and in the numerical simulations of Siviglia et al. (Reference Siviglia, Stecca, Vanzo, Zolezzi, Toro and Tubino2013), as long as the angle stays relatively small.

Figure 9. System of reference used to analyse a geometrically symmetrical bifurcation with non-vanishing angle ${\it\delta}$ .

We consider the symmetrical configuration sketched in figure 1, where downstream anabranches form the same angle ${\it\delta}$ with respect to the main channel. Typical values of ${\it\delta}$ observed in natural gravel-bed bifurcations and in laboratory scale models of braided networks range from $15^{\circ }$ to $30^{\circ }$ (Federici & Paola Reference Federici and Paola2003; Bertoldi & Tubino Reference Bertoldi and Tubino2005; Burge Reference Burge2006; Bertoldi Reference Bertoldi2012).

When the angle is relatively small, the forced solution can be expanded in the following form:

and the general linear solution can be simply obtained by adding the $O({\it\delta})$ component to the free, $O({\it\epsilon})$ , component derived in § 3.

The linear $O({\it\delta})$ problem is solved by considering a rotated coordinate system in channels  $B$ and $C$ , as illustrated in figure 9. Because of the symmetrical geometry of the bifurcation, we can seek a solution for the upper half of the domain and then simply mirror it with respect to the symmetry axis  $y^{A}=0$ . In the lower half the solution for ${\it\eta}$ , $U$ and $D$ is identical, while $V$ is opposite in sign.

In the rotated coordinate system the solution of the order $O({\it\delta})$ for each channel is formally identical to the free solution derived in § 3.2. However, matching the three solutions at the node no longer implies a simple correspondence of the individual components of the velocity vector, $U$ and $V$ , due to the different orientation of the channel axis in the two anabranches. Projecting the velocity components of channel  $B$ in the Cartesian system of channel  $A$ gives the following matching conditions at $x=0$ :

which reduce to those adopted in the analysis of the free response as ${\it\delta}$ vanishes.

Substituting (4.12) into (4.13) we obtain, at $O({\it\delta})$ ,

The second term on the right-hand side of (4.14b ) embodies the forcing effect induced by the deviation of the longitudinal flow direction due to the angle ${\it\delta}$ . It represents the main novelty with respect to the unforced problem and makes the linear system of matching conditions at the node no longer homogeneous. The reader can readily prove that such an additional term generates a component of the overall solution, say in the channel  $A$ , which is odd for the transverse velocity $V$ and even for ${\it\eta}$ , $U$ and $D$ .

Consequently, the forcing effect due to a symmetrical angle between the distributaries produces a symmetrical additional component for bed topography which is unable to affect the bifurcation balance and therefore does not modify the criterion for balanced/unbalanced bifurcations set in terms of the resonant width ratio. This is consistent with the results reported in figure 7, which show that the discharge ratio in the distributaries depends on ${\it\beta}_{A}$ being larger or smaller than ${\it\beta}_{R}$ , irrespective of the presence of a bifurcation angle. Similarly, any symmetric perturbation that does not advantage one of the two branches, such as an even variation of the slope of the distributaries or an equal variation of the downstream boundary condition (in the case of finite length distributaries) and also does not affect the equal partition of the flow in sub-resonant conditions or the unbalance which arises when ${\it\beta}_{A}>{\it\beta}_{R}$ .

The $O({\it\delta})$ solution can be readily obtained on substituting (3.27), (3.19c,d ) and (4.4) into (4.14). Expanding the forcing term in Fourier series we can set the matching conditions at the node in the following form:

where the 1-D component is omitted, as it balances at the node, and $B_{n}$ are the coefficients of the Fourier sine expansion of the second term on the right-hand side of (4.14b ).

The reader will note that the conditions (4.15a,b ) only involve even transverse modes with $m>1$ , so that its structure is valid under both sub- and super-resonant regimes (see figure 4). The solution is straightforward because the different Fourier modes do not interact each other and therefore they can be separately computed. For each even mode $m$ four unknowns $\tilde{{\it\eta}}_{m4}^{A},\,\tilde{{\it\eta}}_{m1}^{B},\,\tilde{{\it\eta}}_{m2}^{B},\,\tilde{{\it\eta}}_{m3}^{B}$ need to be calculated and four matching conditions are available:

As a further simplification, we note that the coefficients $B_{n}$ are given by (4.14b ):

and therefore only the $m$ modes $2,6,10,\ldots$ are actually non-vanishing.

The analytical solution has been shown to provide a consistent agreement with the numerical solution proposed by Edmonds & Slingerland (Reference Edmonds and Slingerland2008) that refers to sub-resonant conditions (Redolfi Reference Redolfi2014, pp. 251–252).