Defining the length parameter in river bifurcation models: a theoretical approach

One‐dimensional models for river bifurcations rely on a nodal point relation that determines the distribution of sediments between the downstream branches. The most widely‐adopted nodal point relation describes the two‐dimensional topographic effects exerted by the bifurcation by introducing two computational cells, located just upstream the bifurcation node, that laterally exchange water and sediments. The results of this approach strongly depend on a dimensionless parameter that represents the ratio between the cell length and the main channel width, whose value needs to be empirically estimated. Previous works proposed calibrating this parameter on the basis of more complete two‐dimensional linear models, which directly solve the momentum and mass conservation equations. This study demonstrates that a full consistency between the one‐dimensional approach and the two‐dimensional models can be directly achieved by adopting different scaling for the bifurcation cell length, which results in a theoretically‐defined and constant dimensionless length parameter. Comparison with experimental observations reveals that this physically‐based scaling yields more accurate predictions of bifurcation stability and discharge asymmetry. This indicates that the proposed method may provide a more reliable and precise estimation of the cell length, potentially improving the performance of one‐dimensional models for bifurcation processes in rivers.


| INTRODUCTION
River bifurcations are found in many natural river settings, such as braiding and anastomosing rivers, alluvial fans, meandering rivers and deltas (Ashmore, 2013;Kästner & Hoitink, 2019;Kleinhans et al., 2013;Marshall & Wohl, 2023). Moreover, they often appear as the result of engineering interventions such as local widening, re-opening of side channels, and construction of longitudinal training dams (de Ruijsscher et al., 2020;Le et al., 2018). Due to their capacity to control the distribution of water and sediment between the downstream branches, bifurcations play a crucial role in determining the morphological trajectory of multi-thread rivers.
One-dimensional modelling approaches have been widely employed to investigate the dynamics of bifurcations in both riverine and deltaic environments (Edmonds et al., 2022;Slingerland & Smith, 2004). The key ingredient of 1-D morphodynamic models is the so-called nodal point relation, which synthetically describes how sediments are diverted at the bifurcation, depending on flow and channel conditions. For example, the simple empirical relation proposed by Wang et al. (1995) enabled for capturing the different morphological trajectories of river bifurcations, highlighting how they can either maintain a balanced discharge partition or progressively lead to the abandonment of one of the branches. However, a more complete modelling of the physical processes is offered by the widely-used relation proposed by Bolla Pittaluga et al. (2003). This relation describes the 2-D morphodynamic effects produced by the bifurcation, through the consideration of two computational cells, located just upstream the bifurcation node, that can laterally exchange water and sediment fluxes. Specifically, the exchange of sediment between the cells is assumed to depend on both the transverse water flux and the lateral bed elevation gradient, which causes sediment deflection through a mechanism of gravitational pull.
A key, general outcome of these models is that the bifurcation evolution is primarily controlled by the channel aspect ratio β (e.g., Bolla Pittaluga et al., 2015;Redolfi et al., 2019), which is conventionally defined as half the width-to-depth ratio of the main channel.
When this parameter exceeds a critical value β C , the initially symmetric and balanced solution becomes unstable, leading to a transition towards a new equilibrium state where a dominating branch receives the majority of the water and sediment fluxes, as illustrated in Figure 1. The critical aspect ratio highly depends on the length of the cells, which is parametrized as L ¼ αW a , where W a is a measure of the main channel width. The dimensionless cell length parameter α is known to be of order one, but its precise value is unknown. Therefore, an empirical calibration is necessary.
To overcome these limitations it is possible to rely on results from fully 2-D models. Specifically, the linear analysis proposed by Redolfi et al. (2016) predicts that instability of balanced configurations occurs when channel aspect ratio exceeds a so-called resonant threshold β R (Blondeaux & Seminara, 1985), which can be calculated by perturbing shallow water and Exner equations. This approach provides predictions of marginal stability conditions that significantly differ from those resulting from the BRT model. This discrepancy is illustrated in Figure 2, from which we can see that for both gravel and sand bed cases, values of β C (left panels) and β R (right panels) are different and show a different dependence on parameters. Notably, resonant aspect ratio predicted by 2-D analysis highly depends on variations of the Chézy coefficient, while the critical aspect ratio given by the BRT approach is nearly unaffected.
The two-dimensional approach can be considered as physicallybased, in the sense that results are derived from the equations of Newtonian dynamics through a series of well-defined and testable assumptions. Moreover, it incorporates a more complete description of the physical processes, and it does not require the calibration of any specific parameter. As a result, it is expected to provide more accurate predictions of marginal stability conditions, as confirmed by Redolfi et al. (2019). This suggests that this approach is generally preferable. However, besides being intrinsically more complicated, the 2-D analysis is not completely alternative to BRT, because the linearized character of the solution prevents from the possibility to study nonlinear phenomena such as interactions between free and forcing effects, and transitions between different equilibrium states.
Therefore, the problem that arises is how to make results of the BRT approach consistent with the physically-based 2-D analysis, preserving the advantages on the 1-D approach in terms of simplicity and ability to explore nonlinear phenomena. The most apparent solution to this problem is to calibrate the bifurcation length parameter to obtain the same critical aspect ratio, as suggested by Redolfi et al. (2019). This calibration process resulted in significantly different values of α depending on the flow parameters and closure formulas, as illustrated in Figure 7 of Redolfi et al. (2019). However, this procedure appears awkward and the final result lacks a direct physical F I G U R E 1 (a) Sketch of the model proposed by Bolla Pittaluga et al. (2003), where L is the bifurcation cell length, while W i and D i indicate the width and the water depth of the individual channels i ¼ fa, b, cg. (b) Variation of the discharge asymmetry ΔQ ¼ ðQ b À Q c Þ=Q a with the channel aspect ratio β ¼ W a =ð2D a Þ, showing the emergence of a "pitchfork bifurcation" of the equilibrium diagram at the critical point β ¼ β C . Please notice that the channel aspect ratio is conventionally defined as half the width to depth ratio of the main channel a. [Color figure can be viewed at wileyonlinelibrary.com] interpretation, to the point that subsequent works have continued to employ an arbitrarily constant value of the parameter α.
In this manuscript, we propose a more straightforward and elegant approach to achieve full consistency between the BRT approach and 2-D linear analysis, based on properly scaling the bifurcation cell length. Moreover, we analyse the implications of this approach for predicting the response of bifurcation stability and degree of asymmetry under varying conditions.

| PHYSICALLY-BASED ESTIMATION OF THE CELL LENGTH
In the original formulation by Bolla Pittaluga et al. (2003), the length of the cells is considered to be proportional to the channel width, namely: F I G U R E 2 Comparison between the critical value of the channel aspect ratio according to the BRT approach (β C , left panels) and the resonant aspect ratio resulting from 2-D linear analysis by Redolfi et al. (2016)  where the coefficient α is considered to be of order one. This choice is motivated by the well-established result that channel width essentially controls the longitudinal scale of fluvial sedimentary patterns such as alternate bars and meandering river bends (Seminara, 2010).
Based on laboratory observations on the length of steady bars forming just upstream the bifurcation, BRT also suggested that values of about 2-3 appear appropriate. This is presumably the main reason for which values around 3 have been frequently adopted in the literature (Table 1).
To revisit this approach, we start from discussing the foundations of the BRT method. The adoption of a two-cell schematization is needed to account for the two-dimensional effects occurring just upstream the bifurcation node. From a physical viewpoint, these effects are associated with the formation of steady bars, which produce a lateral deflection of the bedload flux due to a mechanism of gravitational pull (Redolfi et al., 2019). The model considers the conservation of water and sediment mass within the two cells, and incorporates a physically-based relation for the effect of lateral deflection. However, unlike complete models based on the solution the Saint-Venant equations (e.g., Kleinhans et al., 2008), the BRT scheme does not account for the conservation of momentum. It seems therefore evident that the need to specify the empirical parameter α is the price to pay for ignoring the momentum equations.
Following this idea, we assume that the length of the cells has to be comparable with the characteristic length scale arising from the momentum balance, i.e. the length scale at which inertial effects play a significant role (e.g., Canestrelli et al., 2014). This scale can be quantified by comparing the intensity of flow acceleration and bottom friction terms that appear in the depth-averaged shallow-water equations: where D is the water depth and c is the dimensionless Chézy coefficient, i.e. the ratio between the depth-averaged flow velocity U and the friction velocity u * . Denoting with˜the order of magnitude of the variables, the length scaleL at which the two terms become comparable reads:DŨŨ The resulting parameter c 2 D, which has been already recognized as a significant control for the spacing of fluvial forms (e.g., Ragno et al., 2022;Struiksma et al., 1985), provides a characteristic length that can serve as physically-founded scale for the bifurcation cell length. Following this idea, we introduce a dimensionless coefficient At a first sight, the choice between the scaling given by Equations (1) or (4) seems totally interchangeable, as the associated dimensionless parameters are clearly related as follows: However, here we argue that fixing the dimensionless parameter α 0 is more appropriate, as this allows for directly reconciling results from 1-D and 2-D models. This can be seen by comparing the expressions for the critical and the resonant aspect ratio. According to the BRT model the critical aspect ratio is given by: where Φ D and Φ T are coefficients that take into account the sensitivity of the sediment transport rate to variations of water depth and Shields parameter, c D measures the sensitivity of the Chézy coefficient to variations of water depth (see Redolfi et al., 2019), and r is a dimensionless parameter needed to quantify the effect of lateral slope on bedload transport (Baar et al., 2018;Ikeda, 1982). In this paper, we consider a value r ¼ 0:5, as often assumed in the literature (see Table 1).
The simplified expression for the resonant aspect ratio provided by Camporeale et al. (2007) shows a similar structure: T A B L E 1 Values of the dimensionless cell length parameter α and adopted in the literature for different applications and developments of the BRT approach, and associated values of the dimensionless parameter r that quantifies the effect of lateral bed slope on bedload transport (Baar et al., 2018;Ikeda, 1982). with the noteworthy difference that all the terms on the right hand side are elevated to the power 1/2. To obtain a full agreement between the two models, critical and resonant conditions should clearly match, as they both represent the point at which a pitchfork bifurcation of the equilibrium diagram occurs (Figure 1b). This is ensured by imposing β C ¼ β R , which gives the following value for the bifurcation length parameter: which results in values approximately ranging from 2 to 12, depending on Shields stress and roughness coefficient (see Figure 7 of Redolfi et al., 2019). Noteworthy, this variability is eliminated when computing the re-scaled coefficient α 0 by combining Equations (8) and (5), which simply gives: Therefore, the dimensionless parameter α 0 turns out to be constant and theoretically-defined, at least in near-critical conditions. By extension, it seems reasonable to assume the same value of α 0 to be valid across the entire space of parameters, including the super-critical regions where stable unbalanced configurations are expected to develop. This point will be discussed below in Section 1.
It is worth noting that the above analysis does not explicitly account for the effect of stage-related changes in dune geometry on the Chézy coefficient (Jing & Qin, 2023). Therefore, resolving the cases illustrated in Figure 2c requires a slightly different procedure, which is explained in detail in Appendix A.
In the following, we will examine how this result translates into a distinct response of the bifurcation to varying parameters, analysing the implications in terms of bifurcation stability and degree of asymmetry.

| Bifurcation stability
The first implication of fixing a theoretical value for the scaled coefficient α 0 is that it aligns the critical aspect ratio given by the BRT approach with the resonant value. This eliminates both the qualitative and the quantitative differences in the response to varying parameters highlighted in Figure 2, in the sense that the β C illustrated on the left panels would be exactly equal to the β R appearing on the right.
Since the resonant threshold has been found to provide a satisfactory criteria for assessing the stability of bifurcations (Redolfi et al., 2016), we expect that our scaling approach can enhance the predictive capacity of the BRT model. To test this predictive power, here we analyse the effect of different choices of the length parameter on the F I G U R E 3 Comparison between asymmetry observed in laboratory (Bertoldi & Tubino, 2007), numerical (Siviglia et al., 2013) and field (Zolezzi et al., 2006) bifurcations and theoretical predictions obtained by considering fixed values α ¼ 3 (a), α ¼ 4:5 (b) and α ¼ 6 (c), and fixed α 0 ¼ π 2 =16 (d). Filled markers indicate bifurcations that are observed to be balanced, which should theoretically fall in the sub-critical region, i.e. below the dashed line that indicates the marginal stability condition β ¼ β C . ACC and BA denote Accuracy and Balanced Accuracy metrics, respectively.
From model outcomes we anticipate that under sub-critical conditions (i.e. β < β C ) the observed discharge partition should be balanced, while in the super-critical regime we expect to observe unbalanced configurations. To test this prediction, we first consider classic indicators of classification performance, namely the accuracy ACC ð Þand the balanced accuracy BA ð Þ(see Tharwat, 2018), which are defined as follows: where N indicates the number of unbalanced and super-critical (N SUP U ), balanced and sub-critical (N SUB B ), balanced and super-critical (N SUP B ), unbalanced and sub-critical (N SUB U ) cases, whose sum equals the total number of cases N TOT . Figure 3 show that, when considering the classic scaling, a value near α ¼ 6 yields the best predictions, consistently with the original estimation provided by Bertoldi & Tubino (2007). Nevertheless, an even better prediction is provided by our definition of the cell length parameter α 0 , with the key difference that this method does not need any specific calibration.

Results reported in
It is also worth noting that the same result can be also obtained by choosing α ¼ α R (Redolfi et al., 2019), which according to Equation (8) also ensures β C ¼ β R . However, this option cannot be considered equivalent, for two main reasons. Firstly, the value α R is not constant but needs to be calculated for each experimental point. Secondly, fixing α ¼ α R gives a notably different response of the equilibrium asymmetry when the channel aspect ratio exceeds the critical threshold, as will be highlighted below.

| Bifurcation asymmetry
A direct consequence of fixing the dimensionless parameter α 0 is the change in the response of the bifurcation to increasing values of the channel aspect ratio. Specifically, in an ideal experiment where channel width is increased while keeping the other hydraulic parameters (channel slope, water depth, Shields number) unchanged, the classic scaling (1) implies that L increases proportionally with the channel width, while the proposed scaling gives a constant cell length. From a different viewpoint, fixing α 0 implies that the equivalent value of α decreases as the channel aspect ratio β increases, as indicated by Equation (5). This causes a reduction of the stabilizing effect of the gravitational pull on bedload transport, which explains the sharper response of the bifurcation asymmetry to variations of channel aspect ratio that is illustrated in Figure 4. A direct consequence is that the condition for which one the bifurcation is so asymmetric that the shallower channel is no longer able to transport sediment becomes easier to reach. As a result, unbalanced solutions with both branches active reduce to a narrower range of β values.
The ability of the model to reproduce the observed discharge asymmetry is then tested using laboratory data from Bertoldi & Tubino (2007). The comparison between model predictions and measurements illustrated in Figure 5 reveals that adopting the classic scaling leads to modest variations of the predicted discharge asymmetry with respect to the measured ΔQ. As a result, the value α ¼ 6 that better captures marginal stability conditions (Figure 3), leads to an underestimation of the observed asymmetry, especially for the largest values of the measured ΔQ. The same occurs when employing the value α ¼ α R (different for each experiment) that ensures β C ¼ β R . Conversely, fixing α 0 ¼ π 2 =16 provides a more centered prediction of the observed ΔQ, as highlighted by the lower values of RMSE. This is in line with the analysis conducted by Bertoldi & Tubino (2007), who suggested that a variable value for α should be employed. Specifically, they found that for unstable bifurcations the optimal value of α should decrease with β, which is consistent with our Equation (5).

| DISCUSSION
In the previous Section 1, we have seen how the proposed scaling makes results from the BRT approach fully consistent with the 2-D F I G U R E 4 Equilibrium solutions for a free bifurcation (a) and for a bifurcation forced by a slope asymmetry ΔS ¼ ðS b À S c Þ=S a ¼ 0:01 (b). Gray lines: classic scaling (fixed α parameter); blue lines: new proposed scaling (fixed α 0 parameter). Solid and dashed lines are used to indicate stable and unstable equilibrium solutions, respectively. Adopting the reduced parameter ðβ À β C Þ=β C makes the plot fully independent of the choice of the actual values of α and α 0 . [Color figure can be viewed at wileyonlinelibrary.com] analysis by Redolfi et al. (2016), providing an explanation for the qualitative and quantitative differences in the response of the critical aspect ratio that are illustrated in Figure 2.
In this way, the empirical estimation of the cell length parameter is no longer required. From a practical point of view, this allows for a more accurate estimation of bifurcation stability and equilibrium discharge asymmetry, at least in the explored range of conditions. More fundamentally, our analysis suggests the need of adopting a different scaling length, proportional to the water depth. This has significant implications in the identification of the key dimensionless controls.
For instance Ragno et al. (2021) highlighted that in a channel loop where a bifurcation interacts with a downstream confluence, the effect of the bifurcation is controlled by the parameter R, which is expressed as: This shows how the proposed scaling implies an inverse quadratic (rather than inversely proportional) dependence on the channel aspect ratio β, which is characteristic of alternate bar patterns (Redolfi, 2021).
In general, the scaling length c 2 D brings a precise physical meaning, as it quantifies the length at which flow inertia and friction becomes comparable (e.g., Canestrelli et al., 2014). For this reason, it is not exclusive of bifurcation processes, as it has been recognized to be an important parameter in determining the spacing of various fluvial and deltaic patterns, including bar-pools units, meander bends and bifurcation-confluence loops (Camporeale et al., 2005;Fagherazzi et al., 2015;Johannesson & Parker, 1989;Ragno et al., 2022;Struiksma et al., 1985).
In the following, we will discuss the possible variation of the cell length parameter when the channel aspect ratio departs from the critical value. Finally, we will address limitations and potential future developments.

| Application to non-critical conditions
The constant value for the dimensionless coefficient α 0 is obtained from the analysis of near-critical states, where the system is close to marginal (in)stability conditions. The assumption of extending the applicability of this theoretically-derived value to the entire space of parameters seems quite natural, and is substantiated by the comparison with experimental data presented in Figure 5. However, a theoretical reasoning, based on the theory of morphodynamic influence originally proposed by Zolezzi & Seminara (2001), can provide additional support for this assumption.
Specifically, it is worth considering that, as highlighted by Redolfi et al. (2016), any two-dimensional effect exerted by the bifurcation exponentially decays in the upstream direction. Specifically, at a F I G U R E 5 Comparison between the equilibrium discharge asymmetry ΔQ measured in the laboratory experiments by Bertoldi & Tubino (2007) and values predicted by the BRT model, the latter obtained by considering: fixed α (panels a, b), α ¼ α R (c), and fixed α 0 ¼ π 2 =16 (d). RMSE denotes the root mean square error of the model's predictions. [Color figure can be viewed at wileyonlinelibrary. com] F I G U R E 6 Scaled cell length L=W a as a function of distance from marginal stability conditions, depending on whether α of α 0 are fixed. The parameter L d =W a indicates the (scaled) distance at which the upstream influence exerted by the bifurcations decays by a fraction d, which is expected to provide a upper limit for the bifurcation cell length. Parameters are θ a ¼ 0:06, c a ¼ 11:2. [Color figure can be viewed at wileyonlinelibrary.com] distance L d from the bifurcation, any bed disturbance is expected to decay by a fraction d that is given by: where λ is a dimensionless damping rate, which depends on the parameters β, θ a , c a . At resonant conditions, the damping rate is zero, while for larger values of the channel aspect ratio the damping rate increases (Zolezzi & Seminara, 2001). As a result, the scaled decay length L d =W a depends on the distance from marginal stability conditions, as illustrated in Figure 6 by considering different values of the decay d.
Considering that L d determines the distance at which the upstream influence of the bifurcation can propagate, it should provide an upper bound for the cell length L. This idea originates from Bertoldi & Tubino (2007), who noticed that cell the length parameter L=W a needed to fit the experimentally observed discharge asymmetry decreases when β > β R , and is supported by Redolfi et al. (2016), who highlighted how this calibrated length parameter matches the trend of L d =W a . However, despite its theoretical relevance, the idea that bifurcation cell length should be limited to L d appears practically unfeasible, mainly due to the absence of criteria to select a threshold value for the decay d.
Nonetheless, from the perspective of this work, the important point is that when fixing α 0 the need to limit the cell length becomes less significant than when considering a fixed value of α. This is because the α 0 scaling implies that L=W a automatically reduces when increasing ðβ À β R Þ=β R , as illustrated in Figure 6. Such decrease is ultimately responsible of the sharper response illustrated in Figure 4, and it determines the improved fitting of experimentally-observed discharge asymmetry that is illustrated in Figure 5. Therefore, assuming a constant value of α 0 across the entire parameters space appears practically valid and theoretically more accurate the conventional, fixed α approach.

| Open questions and future perspectives
The above framework appears to be theoretically solid and consistent with existing laboratory and numerical experiments. However, there are still important questions that require further investigation. One of these questions pertains to the application of this approach when the width of the downstream channels is significantly different, as often the case in natural river bifurcations (e.g., Hariharan et al., 2022). In this case, defining the geometry of the bifurcation cells becomes nontrivial, and dedicated experiments (e.g., Islam et al., 2006) should be Another relevant question concerns the role of suspended sediment transport. Since suspended transport is weakly affected by gravitational pull, its effect can be incorporated into the BRT scheme by considering the contributions of the suspended and the bedload fraction separately, as proposed by Iwantoro et al. (2021). However, other factors such as channel width variations, different bifurcation angles, secondary flow, sediment sorting, and influence of tides may also play a role in determining the stability of bifurcations (Kästner & Hoitink, 2019;Marra et al., 2014;Sassi et al., 2013;Slingerland & Smith, 1998). In this sense, further analysis is needed to understand whether these effects can be adequately described within a simple, two cells scheme or whether a more detailed modelling of the flow field, taking into account the actual bifurcation geometry, is required (Kästner & Hoitink, 2020).
Finally, it is worth mentioning that the proposed scaling ensures that the critical aspect ratio matches the approximated value of the resonant aspect ratio, obtained from the expression by Camporeale et al. (2007), which is based on neglecting water surface effects. However, this seems a minor limitation, considering that the error of the simplified expression is of the order of a few percent (see Appendix B).

| CONCLUSIONS
In this study, we have addressed the problem of estimating the bifurcation length parameter that arises in one-dimensional bifurcation models based on the widely-used, two-cell approach proposed by Bolla Pittaluga et al. (2003). The analysis has demonstrated that by considering a specific scaling for the cell length parameter, full consistency between 1-D approaches and 2-D linear analysis can be achieved. Specifically, the scaled length parameter α 0 , defined as the ratio between the cell length L and the quantity c 2 a D a , is found to assume a constant, theoretically-defined value α 0 ¼ π 2 =16. Therefore, the cell length to be implemented in one-dimensional bifurcation models is simply given by the following expression: or, equivalently, by setting the classic parameter α of BRT-based nodal point relations as follows: where β is half the width to depth ratio, D a is the water depth in the main channel, and c a is the dimensionless Chézy coefficient (i.e. the ratio between the average flow velocity and the friction velocity).
The results of this study indicate that adopting this scaling leads to a quantitatively and qualitatively different response of the bifurcation to changing parameters. In particular, a key difference is the stronger effect of the channel width-to-depth ratio, resulting in a sharper variation of the discharge asymmetry when the critical threshold condition is exceeded.
Furthermore, the analysis of existing data from numerical and laboratory experiments has revealed that the proposed approach provides a better prediction of both bifurcation stability and the asymmetry of equilibrium configurations observed in laboratory and numerical models of gravel-bed river bifurcations. Therefore, the present method offers reasonable, physically-based estimation of the bifurcation length parameter, particularly in situations where direct empirical calibration is not feasible.

AUTHOR CONTRIBUTIONS
Single-autor paper.

ACKNOWLEDGEMENTS
This study was partially carried out within the PNRR research activities of the consortium iNEST (Interconnected North-Est Innovation Ecosystem) funded by the European Union Next-GenerationEU thus accounting for the form drag exerted by dunes. In this case, the Chézy coefficient can be derived from the following formulae: where the grain-related water depth is given by: As illustrated in Figure 7 the resonant aspect ratio resulting from the simplified expression (7) turns out to be very similar to the value provided by the complete model, showing a maximum error of about 5% for sand bed cases, and 7% for gravel bed channels, the latter falling below 3% when considering sub-critical (i.e. Fr < 1) conditions only.