3.1. The free response

An examination of the differential problem (2.7) suggests that $\mathcal {R}$ fundamentally controls the dynamics of the system. As emphasized by previous theoretical studies (Ragno et al. Reference Ragno, Redolfi and Tubino2021; Ragno, Tambroni & Bolla Pittaluga Reference Ragno, Tambroni and Bolla Pittaluga2022b), the bifurcation parameter, which is inversely proportional to $\mathcal {B}_{0}$, is found to control the number and relative stability of equilibria: if $\mathcal {R}$ is higher than the critical threshold $\mathcal {R}_{cr}$, then flow and sediment fluxes are distributed equally at the bifurcation. In this case, $\Delta D =0$. The situation changes radically when $\mathcal {R}$ is lower than $\mathcal {R}_{cr}$, since the balanced solution becomes unstable, with the system attaining a stable state where one of the two branches captures most of the discharge, i.e. $\Delta D \ne 0$. According to the BRT model, the critical bifurcation parameter is given by

(3.1)\begin{equation} \mathcal{R}_{cr} = \tfrac{1}{4} \left(\varPhi_{T} + \varPhi_{D} - c_{D} - \tfrac{3}{2} \right) , \end{equation}

where coefficients $\varPhi _{T}$, $\varPhi _{D}$ and $c_{D}$, which measure the sensitivity of the sediment transport rate and Chézy coefficient to variations of flow depth and Shields parameter, are defined as

(3.2a–c)\begin{equation} \varPhi_{T}:=\left.\frac{\tau_{*0}}{\varPhi_{0}}\,\frac{\partial \varPhi}{\partial \tau_{*}}\right|_{\tau_{*0}}, \quad \varPhi_{D}:=\left.\frac{D_{0}}{\varPhi_{0}}\,\frac{\partial \varPhi}{\partial D}\right|_{D_{0}}, \quad c_{D}:=\left.\frac{D_{0}}{c_{0}}\,\frac{\partial c}{\partial D}\right|_{D_{0}} , \end{equation}

and their explicit expression depends on the friction and transport formula adopted (Redolfi, Zolezzi & Tubino Reference Redolfi, Zolezzi and Tubino2019).

Taking the viewpoint of condensed matter physics, here it is proposed that the critical bifurcation parameter can be conceived as a morphodynamic analogue of the Curie temperature above which spontaneous magnetization disappears in a ferromagnetic material (e.g. Kadanoff et al. Reference Kadanoff, Gotze, Hamblen, Hecht, Lewis, Palciauskas, Rayl and Swift1967; Wilson Reference Wilson1983). If this analogy is pursued, then it seems reasonable to follow the fundamental idea apparently first introduced by Landau (Reference Landau1937), who suggested to study the behaviour of systems showing a second-order phase transition near a critical point in terms of an order parameter (e.g. Landau & Lifshitz Reference Landau and Lifshitz1980; Domb Reference Domb1996). In this problem, the order parameter arises naturally from system (2.7), coinciding with the depth elevation asymmetry, $\Delta D$.

In order to measure the deviation from the critical point, the parameter $\epsilon$ is introduced:

(3.3)\begin{equation} \epsilon := \frac{\mathcal{R}_{cr}-\mathcal{R}}{\mathcal{R}_{cr}} \ll 1. \end{equation}

From this, several properties of the order parameter follow (Kadanoff et al. Reference Kadanoff, Gotze, Hamblen, Hecht, Lewis, Palciauskas, Rayl and Swift1967). First, $\Delta D = 0$ when $\epsilon < 0$ (subcritical phase), but $\Delta D$ is different from zero when $\epsilon >0$. Second, a discontinuous jump of $\Delta D$ is absent when $\mathcal {R} \rightarrow \mathcal {R}_{cr}$. Third, when $\epsilon > 0$ (supercritical phase), $\Delta D$ can assume two or more values under the same physical conditions. Indeed, either the left or the right branch can capture the largest fraction of water and sediment.

The phase transition, i.e. the evolution of a river bifurcation from a balanced to an unbalanced flux partitioning, can be solved once the functional relation $f(\epsilon,\Delta D)$ is computed. Near criticality, where the order parameter $\Delta D$ is small, each variable, say generically $X$, is assumed to be a slowly varying function in time and analytic in $\epsilon$. Thus the slow time scale $T=\epsilon t$ is defined, and the following expansion of the solution holds:

(3.4)\begin{equation} X = X_{0} + \sum_{j=1}^{3} X_{j} (\epsilon^{\beta})^{j} + O(\epsilon^{4\beta}), \end{equation}

where the value of the expansion coefficient $\beta$ is a priori unknown.

A strategy to determine the order of the critical parameter $\epsilon$ is to follow a heuristic argument adopted widely in the field of hydrodynamic and morphodynamic stability (e.g. Stuart Reference Stuart1960; Colombini, Seminara & Tubino Reference Colombini, Seminara and Tubino1987): reason of symmetry suggests that the linear component $X_{1}$ is modified at the third order, where the term $\epsilon \,{\rm d}(\epsilon ^{\beta } X )/{\rm d}T$ must be balanced by linear and nonlinear terms proportional to $\epsilon ^{3\beta }$. Consequently, the critical index $\beta$ is equal to $1/2$, which corroborates the analogy with Landau's theory, as discussed below.

Substituting (3.4) into system (2.7), and equating like powers of $\epsilon$, a cascade of algebraic problems is obtained. At the order $\epsilon ^{3/2}$, the following nonlinear ordinary differential equation governing the evolution of the order parameter $\Delta D$ is found:

(3.5)\begin{equation} \frac{{\rm d} \Delta D}{{\rm d}T} = \epsilon \alpha_{1}\,\Delta D + \alpha_{2}(\Delta D^{3}), \end{equation}

where the coefficients $\alpha _{1}$ and $\alpha _{2}$ depend on $\epsilon$, and on the specific reference conditions (i.e. $\tau _{*0}$ and $c_{0}$). Their expression is reported in Appendix A along with the computations to derive (3.5).

Equation (3.5) is of the Stuart–Landau type, and admits an analytical solution of the form

(3.6)\begin{equation} \Delta D(t)= \left[ - \frac{\alpha_{2}}{\alpha_{1}} + \left(\Delta D(0)^{-2} + \frac{\alpha_{2}}{\alpha_{1}} \right) \exp(- \epsilon \alpha_{1} t) \right]^{-1/2} . \end{equation}

If the cubic term in (3.5) is neglected, then the order parameter $\Delta D$ exhibits an exponential behaviour, with the growth rate governed by the coefficient $\alpha _{1}$. For $\epsilon >0$, the nonlinear term inhibits the unbounded growth since the algebraic coefficients $\alpha _{1}$ and $\alpha _{2}$ have opposite sign (figure 9). Hence in the first stages, the growth of asymmetry development is governed by the linear component of the solution. In this phase, the solution follows a nearly exponential behaviour that is proportional to $\exp (\alpha _{1} t)$. Then higher-order terms become important, and the curve changes concavity until it reaches the final equilibrium state (figure 4).

Figure 4.

Evolution of the order parameter $\Delta D$ as scaled by its equilibrium value $\Delta D_{eq}$ in a reference supercritical phase. The dashed line denotes the linear solution. Without loss of generality, the plot refers to a gravel-bed river for which the Parker (Reference Parker1978) transport formula is used.

Thus as $t \rightarrow \infty$, an equilibrium solution $\Delta D_{eq}$ is reached asymptotically:

(3.7)\begin{equation} \Delta D_{eq} ={\pm} \sqrt{- \frac{\alpha_{1}}{\alpha_{2}}}\,\epsilon^{1/2} . \end{equation}

The equilibrium solution (3.7) describes mathematically a supercritical pitchfork bifurcation (figure 5a). When $\epsilon <0$, the only real solution is $\Delta D_{eq} = 0$, i.e. the balanced bifurcation. Conversely, two real solutions exist when $\epsilon >0$, showing physically the steady state where one of the two branches captures most of the water and sediment fed at the bifurcation.

Figure 5.

(a) Phase diagram of a free bifurcation. The order parameter is plotted as a function of $\epsilon$. The yellow star indicates the critical point. (b) Potential function $\mathcal {G}$ against the steady order parameter $\Delta D_{eq}$ for two values of the scaled bifurcation parameter $\epsilon$. Without loss of generality, both plots refer to a gravel-bed river for which the Parker (Reference Parker1978) transport formula is adopted.

From (3.7), it follows that in the supercritical phase, the order parameter is proportional to the scaled bifurcation parameter $\epsilon ^{\beta }$, with $\beta = 1/2$. Despite the ratio between $\alpha _{1}$ and $\alpha _{2}$ in general not being constant, the critical exponent keeps the same value. This observation highlights an essential property of the system, namely universality. Such universal behaviour is a typical feature of phase transitions (e.g. Kadanoff Reference Kadanoff1990). The critical exponent conforms with the value obtained from the mean-field theory employed by Landau (Reference Landau1937). In order to obtain the latter result, Landau followed a classic rule of thermodynamics whereby the most probable value of any (macroscopic) thermodynamic variable is the one that minimizes the free energy. Following the thread of Landau, a morphodynamic potential $\mathcal {G}=f(\epsilon,\Delta D; \tau _{*0},c_{0})$ can be defined as

(3.8)\begin{equation} \mathcal{G} =- \left[ \tfrac{1}{2}\epsilon\alpha_{1}\,\Delta D^{2} + \tfrac{1}{4}\alpha_{2}\,\Delta D^{4} \right], \end{equation}

where no odd powers of $\Delta D$ are present due to the symmetry requirement. In other words, the morphodynamic potential cannot depend on the sign of $\Delta D$ if no external factors are included.

Figure 5(b) shows the morphodynamic potential of a reference gravel-bed river for which a typical threshold transport function is employed. Overall, the residual effect of Chézy coefficient $c_{0}$ on $\alpha _{1}$ and $\alpha _{2}$ is very weak, to the point that it vanishes exactly when adopting a friction formula of the Manning–Strickler type (e.g. Ferguson Reference Ferguson2007; Parker et al. Reference Parker, Wilcock, Paola, Dietrich and Pitlick2007). Similarly, there is no residual effect of the Shields number $\tau _{*0}$ when considering sand-bed cases, where the transport function is usually estimated through a power-law relation, while a rather important effect remains in gravel-bed cases. Consequently, for sand-bed rivers, the equilibrium value predicted by the model can be approximated by the expression

(3.9)\begin{equation} \Delta D_{eq}^{sand} ={\pm} 0.55 \epsilon^{1/2}, \end{equation}

which, it is worth stressing, is strictly valid for $\epsilon \rightarrow 0$.

Summarizing, in this subsection it has been shown how river bifurcations can be interpreted as critical phenomena showing several analogies with thermodynamic phase transitions. The definition of the potential (3.8) is particularly useful, since it allows us to elucidate characteristics of the morphodynamic system taking advantage of theories developed in the context of condensed matter physics. Specifically, it will be seen below how the morphodynamic potential proves to be convenient to investigate the response of bifurcations to an (applied) external forcing.

3.2. The effect of an external forcing

Up to this point, the free response of the system has been analysed. However, in nature, bifurcations are seldom free, with allogenic forcing factors coexisting with the free response. Physically, a forcing effect can be caused by a slope advantage to one of the two branches (e.g. due to backwater effects, different channel lengths due to meander cut-offs) or by the curvature of the main channel causing a secondary flow component at the node (Redolfi et al. Reference Redolfi, Zolezzi and Tubino2019). In the context of Landau theory, an external yet physically unspecified steady forcing $\mathcal {F}$ is considered here. This is equivalent to considering $\mathcal {F}$ acting on a much slower time scale than the intrinsic time scale of bifurcation evolution (2.6). If this forcing is assumed to be linearly dependent in the order parameter, then the (free) potential (3.8) is modified accordingly as

(3.10)\begin{equation} \mathcal{G} =- \left[\tfrac{1}{2} \epsilon \alpha_{1}\,\Delta D^{2} + \tfrac{1}{4} \alpha_{2}\,\Delta D^{4} - \varphi \mathcal{F}\,\Delta D \right], \end{equation}

which, when minimized, yields

(3.11)\begin{equation} \mathcal{F} = \varphi^{-1} ( \epsilon \alpha_{1}\,\Delta D_{eq} + \alpha_{2}\,\Delta D_{eq}^{3}) , \end{equation}

where $\varphi$ is a calibration coefficient.

One direct consequence of the applied forcing is the non-zero value of $\Delta D_{eq}$ for any given value of $\epsilon$. In other words, there is no spontaneous symmetry breaking of the system (figure 6a). In the subcritical phase, the system is dominated by the forcing effect, with the depth asymmetry varying smoothly with $\mathcal {F}$. In this case, there is just one real stable solution. The picture changes radically in the supercritical phase. As (3.11) indicates, there is a range of values of $\mathcal {F}$ such that three real solutions exist. However, one of these solutions turns out to be unstable.

Figure 6.

The relationship between the free and forced responses. (a) Order parameter as a function of forcing $\mathcal {F}$. (b) Forcing threshold $\mathcal {F}_{t}$ against $\epsilon$. Coloured lines stand for gravel-bed rivers with different values of Shields number $\tau _{*0}$. The dashed line denotes a reference sand-bed river.

The range of values of $\mathcal {F}$ where the solution is no longer single-valued, thus $-\mathcal {F}_{t} < \mathcal {F}< \mathcal {F}_{t}$, can be derived readily by differentiating (3.11) with respect to $\Delta D_{eq}$. The quantity $\mathcal {F}_{t}$ is given by

(3.12)\begin{equation} \mathcal{F}_{t} \sim \left(\alpha_{1}|\epsilon| \right)^{3/2} \alpha_{2}^{-1/2}, \end{equation}

and it provides the threshold value of the forcing whereby the free and forced responses have the same order of magnitude (Landau & Lifshitz Reference Landau and Lifshitz1980). The threshold $\mathcal {F}_{t}$ corresponds to a forcing able to suppress entirely the spontaneous response of the system. According to (3.12), higher values of the threshold are required as long as the deviation of the system from critical conditions increases. The characteristic value of the spontaneous response is provided by (3.7), namely $\Delta D_{eq}^{S} \sim (\alpha _{1}|\epsilon |/\alpha _{2})^{1/2}$, while the forced-induced component of the order parameter is estimated as $\Delta D_{eq}^{F} \sim \chi \mathcal {F}$. The quantity $\chi$ is an equivalent of the magnetic susceptibility, namely a measure of the sensitivity of the order parameter under the action of an external forcing. The morphodynamic susceptibility is described by $\chi \sim (\alpha _{1} \epsilon )^{-1}$, indicating that the system can be highly affected by weak variations of the forcing as long as $\epsilon$ approaches criticality.

The computation of function $\mathcal {F}_{t}$ enables the definition of a morphodynamic bound that divides two qualitatively different behaviours of the bifurcation system when subjected to allogenic factors (figure 6b). When $\mathcal {F} \gg \mathcal {F}_{t}$, the system is forced-dominated, with the free response almost completely suppressed. In this case, the order parameter can be estimated from (3.11), and it is equal to $\Delta D_{eq}^{F} \sim (\mathcal {F}/\alpha _{2})^{1/3}$. On the other hand, when $\mathcal {F} < \mathcal {F}_{t}$, an interaction between the autogenic and allogenic dynamics is promoted, with the possibility of undergoing hysteresis and counterintuitive solutions of the bifurcation system.

Nevertheless, the definition of the forcing might appear abstract to the reader. To overcome this limitation, the recent numerical analysis performed by Redolfi et al. (Reference Redolfi, Zolezzi and Tubino2019), in which different kind of forcing effects are modelled physically, lends itself to being particularly instructive. Assume that the upstream channel is a meander bend. Redolfi et al. (Reference Redolfi, Zolezzi and Tubino2019) have shown that channel curvature, here denoted $\mathcal {C}$, generates a further spiral flow component in the equation of transverse transport (2.2), which produces an uneven distribution of the incoming fluxes at the node. In the present modelling framework, this additional term is absent, since the forcing $\mathcal {C}$ is imposed simply as linearly dependent on the order parameter. However, as shown in figure 7, the analytical solution (3.12) provides a reasonable estimate of the forcing threshold when compared with the numerical solution. It is worth mentioning that in order to make the analytical results consistent with Redolfi et al. (Reference Redolfi, Zolezzi and Tubino2019), a slightly different definition of the critical parameter is adopted. Namely, the deviation from criticality $\hat {\epsilon }$ is defined with respect to the half width-to-depth ratio:

(3.13)\begin{equation} \hat{\epsilon} = \frac{\mathcal{B} - \mathcal{B}_{R}}{\mathcal{B}} , \end{equation}

where $\mathcal {B}_{R}$ is the resonant half width-to-depth ratio, which defines the threshold between prevailing downstream or upstream influence of two-dimensional morphodynamic changes (Zolezzi & Seminara Reference Zolezzi and Seminara2001; Redolfi et al. Reference Redolfi, Zolezzi and Tubino2016). Recalling the physical meaning of morphodynamic resonance, when the channel width-to-depth ratio is smaller (greater) than the resonant value, a regular sequence of damped steady alternate bars is predicted to form downstream (upstream) of a persistent obstacle (Zolezzi & Seminara Reference Zolezzi and Seminara2001). If, as in the present work, the obstacle is a channel bifurcation, then super-resonant conditions correspond to the case of an unbalanced bifurcation, with the resonant width-to-depth ratio related closely to the critical bifurcation parameter discussed previously. Specifically, as suggested by Redolfi et al. (Reference Redolfi, Zolezzi and Tubino2019), the resonant value can be computed by calibrating the bifurcation length $\xi$ to obtain the same critical width-to-depth ratio. It can be realized that the alternative critical parameter (3.13) does not change the qualitative result.

Figure 7.

Comparison between the analytical expression (3.12) of the forcing threshold, which is given for different values of coefficient $\varphi$ (blue lines), and the numerical solution of Redolfi et al. (Reference Redolfi, Zolezzi and Tubino2019) (black line). The forcing is exemplified by the channel curvature $\mathcal {C}$, while $\hat {\epsilon }$ is the scaled half width-to-depth ratio (3.13).