Theoretical Investigation of Equilibrium Dynamics in Braided Gravel Beds for the Preservation of a Sustainable Fluvial Environment

: Gravel bars have an important role in the exchange between surface and subsurface waters, in preventing and mitigating riverbank erosion, in allowing the recreational use of rivers, and in preserving ﬂuvial or riparian habitats for species of ﬁshes, invertebrates, plants, and birds. In many cases, gravel bars constitute an important substrate for the establishment and development of ground ﬂora and woody vegetation and guarantee higher plant diversity. A sustainable management of braided rivers should, therefore, ensure their ecological potential and biodiversity by preserving a suitable braiding structure over time. In the present study, we propose an analytical–numerical model for predicting the evolution of gravel bars in conditions of dynamical equilibrium. The model is based on the combination of sediment balance equation and a regression formula relating dimensionless unit bedload rate and stream power. The results highlight the dependence of the evolving sediment particles’ pattern on the ratio of initial macro-bedforms longitudinal dimension to river width, which determines the gradual transition from advective and highly braiding to diffusive transport regime. Speciﬁcally, the tendency to maintain braiding and ﬂow bifurcation is associated with equilibrium average bed proﬁles and, therefore, equilibrium average stream power characterized by the maximum period that does not exceed transverse channel dimension.


Introduction
Braided rivers are typically organized in multithread patterns, with central or alternate bars dividing channels whose reach-scale statistical properties are controlled by discharge, bed slope, and grain size [1][2][3][4][5]. Braiding in gravel-bed streams may derive from one of several mechanisms, including chute cut-off of point or alternating bars, flow division around a transverse bar, dissection of multiple bars, and central bar formation [6]. Interactions among braided pattern, channel morphology, and bedload transport are significant [7]. Indeed, bedload transport in gravel-bed rivers is a very intricate phenomenon that is characterized by high variability in time and space [8][9][10][11]. Such variability derives from different mechanisms that act at different spatial scales [12][13][14][15], ranging from grain to reach: (1) the movement of single grains (e.g., [16][17][18][19]); (2) the creation of small bed forms, such as sediment waves [6], bed waves [20], and bedload sheets [21,22]; and (3) large-scale morphological processes, such as bar formation and migration, activation and abandonment of anabranches and side-channel shift [6,23,24]. The combination of all these processes results in a fluctuating nature of bedload transport, which was initially recognized by field surveys [25][26][27][28] and later confirmed by innovative measurement procedures [29][30][31][32]. As a matter of fact, several experimental studies performed under steady flow and sediment-supply conditions highlighted huge fluctuations of bedload transport rates both in single-thread [12,13,22,33] and braided channels [6,33]. At grain scale, channel morphological dynamics are the result of single-particle displacements from erosion sites and sediment diffusion in river bedload (e.g., [16][17][18][19]). Specifically, ref. [16] introduced a new conceptual model for the diffusion of moving-bed particles, which suggests that particle motion is diffusive in nature and includes at least three ranges of temporal and spatial scales characterized by different diffusion regimes. Later on, by means of particle tracers, ref. [34] confirmed the diffusive nature of particle displacement in real gravel-bed rivers and investigated gravel particle paths and their contribution to self-forming processes, such as gravel bar development. Particle paths were found to be mostly bank-line parallel, even in river bends. Additionally, transverse bedload was evaluated as 4% of the longitudinal counterpart. Results from physical models of braided rivers (e.g., [24]) highlighted the influence of bars in acting as deposition sites for moving particles, with over 80% of sediment tracers that were found at side, diagonal, and point bars. So far, very few examples of analytical or semi-analytical treatment of sediment continuity equation are available in the literature (e.g., [35]). In any case, none of these provides a solution in terms of actual bed elevation function. A two-dimensional fully numerical model was developed by [36] to simulate bar growth, channel widening, and scour holes at bar lee. In the present study, we make a system of two-dimensional sediment mass balance and a suitable constitutive relationship that relates dimensionless bedload rate and stream power in the case of braided gravel beds at dynamical equilibrium. Such a relationship was derived by the authors in a previous high-grain Reynolds numbers experimental study [3]. The crucial assumption enabling the analytical-numerical solution of the system in terms of actual bed elevation function is represented by the recognition of the critical state of flow along the braids and, as a consequence, by the equivalence of local section-average energy loss and local section-average bed slope (e.g., [3,[37][38][39][40][41]). Ten experiments were carried out by [3] with fluvial gravel (D 50 = 9 mm, D 10 = 5 mm, and D 90 = 13 mm) to analyze the relationship between flow stream power and bedload transport rate in braided rivers at high grain Reynolds numbers. The experimental runs were conducted in a rectangular, 1 m wide, and 16.7 m long flume with Plexiglas walls. The experimental device consisted of a water re-circulating/sediment feeding system. Gravel was added from a 2.3 m long sediment-supplying basin (located upstream of the flume and downstream of water inflow, with imposed average bottom slope) by means of a sediment feeder mounted on a travelling crane. The mixture of sediment and water coming from the supplying basin entered the flume and modified its slope until the achievement of equilibrium conditions. The runs differed from each other in terms of flow rate and initial slopes of flume and sedimentsupplying basin bottom. See Table 1 for a summary of the experimental conditions, where S 0 and S eq indicate initial and equilibrium flume bottom slopes, S 0SB initial supplying basin bottom slope, and t eq the running time needed to achieve equilibrium. Table 1. Summary of the experimental laboratory conditions for the derivation of the equilibrium relationship between dimensionless bedload rate and stream power [3]. The equilibrium grain Reynolds number, which ranged from 601 to 710, was consistent with the condition of fully turbulent flow. Flow rates were measured by an electromagnetic current meter located on the recirculation pipe. For every run, water surface and bed elevations were measured by a manual point gauge mounted on a travelling carriage at 55 flume and 10 supplying basin cross-sections on average. Each of the flume cross-sections, which were spaced 30 cm apart, was monitored over time at 11 grid-points spaced 10 cm apart, with 605 points surveyed at the different sampling times. Sampling times number varied from 3 (R1) to 9 (R9), as a function of the total running time. As a consequence, during each single run, the total number of bed and free surface elevation measurements ranged from 1815 (for R1) to 5445 (for R9). Sediment flux at the flume inlet, which was proportional to the imposed sediment supplying basin bottom slope and water discharge, was estimated as the average volume of sediment eroded by water per unit time from the sediment-supplying basin. The average bedload output rate of the flume was measured by evaluating the cumulative weight of sediment conveyed into a flume downstream tank. The achievement of equilibrium conditions was determined by verifying that input and output bedload rates were equal and constant.
For each run, at equilibrium, section-averaged dimensionless unit bedload rate q * bx (where x indicates stream-wise direction, the over-bar indicates section-average and the asterisk dimensionless quantities) was evaluated at flume scale (which is representative of a fluvial reach) as a function of the section-averaged dimensionless unit stream power ω * [3]. Data were then analyzed and compared to those available from previous laboratory [2,[42][43][44] and field studies [1,[45][46][47]. The analysis highlighted how different combinations of grain Reynolds number and stream power range led to different types of equilibrium correlation laws (linear or power). In particular, for high grain Reynolds numbers (>601) and for ω * < 1, which are the conditions investigated in the present work as typical field conditions, the experimental data could be well-interpreted by the following power law: Equation (1) was obtained by regression from our 10-run experimental laboratory data and field data collected along 66 braided gravel-bed rivers all over the world.
As reported by [3], in the case of braided gravel-bed rivers at high grain Reynolds numbers, Equation (1) estimates bedload transport rates more accurately than other formulas previously appeared in the literature [2,[48][49][50]. Data of water surface elevations exhibited small oscillations about an average value during every run. Free surface oscillations confirmed the establishment of local critical conditions in the braids, as expected from the equilibrium Froude number of the experimental tests, which ranged between 1.15 and 1.33.

Methods
Let function η = η(x, y, t) represent bed elevation above a given reference level, with y and t indicating cross-sectional coordinate and time counted from equilibrium, respectively. We will proceed by subdividing η into section-average (η) and corresponding deviation (η ): In Equation (2) and in what follows, a indicates a constant, i feq the reach-scale equilibrium bed slope, N the de-trended section-averaged bed profile (which reflects large-scale bedforms, such as bars) and η accounts for small-scale bedforms, such as sediment waves or sediment sheets. Additionally, where B indicates total bed width. Unit-width sediment continuity reads: where p is sediment porosity, s indicates a source term (possibly attributable to slope or tributary sediment inputs) and q b the bedload rate vector per unit channel width: with: and q by = q by (x, y, t) ∼ = q by (x, y, t) Note that the hypothesis of negligible average transverse bedload and total transverse bedload practically coinciding with its deviation, expressed by Equation (7), is consistent with the experimental analysis carried out by [34], which highlighted how, in braided rivers, transverse bedload is only a small percentage of the longitudinal counterpart. The substitution of Equation (2) and Equations (5)- (7) into Equation (4) followed by sectionaveraging yields: with: The difference of Equations (4) and (8) is: where: As above summarized, from their previous laboratory experiment [3], the authors had found that in braided gravel beds at high grain Reynolds numbers and river reachscale equilibrium: where the average unit stream power is expressed by: In Equations (12) and (13) and in what follows, d indicates representative grain size, Q flow rate, g acceleration due to gravity, ρ water density (with γ = ρg), ρ s sediment density, and j section-averaged energy loss. It is worth recalling here that the term "reachscale equilibrium" refers to a situation in which the input bedload rate to the reach is equal to the output bedload rate from the reach, and the reach-scale bed trend slope is equal to the reach-scale free surface slope. Since in braided gravel beds water flows in critical conditions within each of the braids (e.g., [3,[37][38][39][40][41], the corresponding sectionaveraged specific energy E s , given by the sum of section-averaged kinetic energy U(x) 2 /2g (with U(x) indicating section-averaged velocity) and section-averaged flow depth H(x), is characterized by zero-mean longitudinal derivative (e.g., [51]): where b indicates the equilibrium-constant total wetted width and Fr section-averaged Froude number (Fr = U/ gH). In these conditions, section-averaged energy loss is equal to section-averaged bed slope: Therefore: and with: Note that values of j < 0 do not constitute a physical contradiction in terms of energy valley-downstream variation because section-averaged bed elevation η does not have a physical but just a mathematical meaning. The same value of η can be representative of completely different transverse η-distributions. Water cannot be on top of bars and only flows within highly meandering braids. Equation (12) may also be thought as deriving from a formally analogous equation relating the actual (total) quantities: In this case, (20) with q bx given by Equation (17) and By analogy, we will assume (and we will justify later) that: Substitution of Equation (17) into Equation (8) and Equations (21) and (22) into Equation (10), respectively, gives a variable-coefficient diffusion equation and a variablecoefficient advection-diffusion equation with source terms: and ∂η ∂t Note that we are consistently using the term diffusion equation and advectiondiffusion equation for (23) and (24), respectively, because both N and η represent a length proportional to the section-averaged or local concentration of sediment particles. A linearization of Equations (23) and (24) will be achieved here by taking the reach-averaged values of ∂N/∂x, ∂η /∂x and ∂η /∂y in the parentheses on their right-hand sides. The reason for doing this is that, given the experimentally detected wave-like behavior of functions N and η , the parentheses on the right-hand side of Equations (23) and (24) represent a sort of oscillatory diffusion coefficients, which can assume positive as well as negative values depending on the sign of local and section-average bed slope. Considering a generalized version of Fick's law, "oscillatory diffusion coefficients" means that at small (in the case of ∂N/∂x) and very small (in the case of ∂η /∂x and ∂η /∂y) scale, there is a pseudo-periodical alternation of zones characterized by a net flux of sediment towards higher concentrations and zones characterized by a net flux of sediment toward lower concentrations, which overall creates a mutual compensation. Note that the same type of compensation does not take place in the case of the oscillatory velocity (24), because advection moves the particles regardless of their local concentration. Based on the experimental data collected during the experimental runs performed by [3] and, as expected, due to the increasing frequency of decreasing-scale components of the bed elevation function, it results that: where the over-bar and subscript R jointly indicate reach-average. As an example, in the case of the reference run R4 (see Table 1), i f eq = 0.015, ∂N/∂x R = −0.002, The conditions expressed by Relationship (25) enabled the following approximation of Equations (23) and (24), respectively: and ∂η ∂t It should be emphasized that the linear (and, in the above perspective, leading) parts of Equations (21) and (22) can respectively be re-written as: and where τ * b indicates dimensionless section-averaged bed shear-stress: and ε r relative bed roughness (ε r = d/H). Note also that Equation (29) is fully compatible with that given in more general cases for mild secondary currents by [52] and later by [53], when put in the form: In Equation (31), α is a dimensionless constant (equal to 4.93 in the study by [53]), and r is a coefficient given by the ratio of lift to drag forces. In Tealdi's et al. model [53], q sy is the cross-sectional bedload rate along the sloping and counter-sloping side-wall (with the slope expressed by ∂η/∂y) and τ * s the dimensionless shear-stress at the middle bank of a trapezoidal cross-section. In the present model, τ * b is obtained as a function of the sectionaveraged flow depth, which is practically measured starting from the average braided-bed level; additionally, q by ∼ = q by represents the (zero-mean) cross-sectional bedload rate along the sloping and counter-sloping braid walls. Therefore, Equations (29) and (31) may be considered as coinciding with α = 0.2 and r, which represents the hydrodynamical engine of bedload transport, consistently given by the ratio of Froude number squared (a measure of lifting potential magnitude) to relative bed roughness squared (a measure of drag potential magnitude).
Equations (26) and (27) are two coupled partial differential equations. The first is a diffusion equation with constant diffusion coefficient D = 2Γi f eq /(1 − p) and space-time dependent source term S(x, t)/(1 − p). The second is an advection-diffusion equation with constant diffusion coefficient Variable-coefficients advection-diffusion equations were already solved by the authors in a stochastic framework in the case of river-flow transport of dissolved chemicals or fine suspended sediments in the presence of stationary bed morphology heterogeneity [54][55][56]. In the present case, due to the intrinsically non-stationary nature of bed morphology heterogeneity in braided rivers (e.g., [3]), no analytically treatable stochastic approach could apply. Note that (26) and (27) jointly reproduce the diffusive nature of particle displacement in gravel beds highlighted by [16] and later experimentally verified by [34] in the case of real gravel-bed rivers. From (26) and (27), and in the absence of external forcing (S(x, t) = 0, s (x, y, t) = 0), one can write: The meaning of Equation (32) is that, besides the isotropic diffusion undergone by the small-scale bed undulations represented by η , the latter experience upstream migration where section-averaged bed elevation tends to decrease (∂N/i f eq ∂t = [2Γ/(1 − p)]∂ 2 N/∂x 2 < 0) and downstream migration where the section-averaged bed elevation tends to increase (∂N/i f eq ∂t = [2Γ/(1 − p)]∂ 2 N/∂x 2 > 0). As an elementary example, let us start from a section-averaged bed elevation given by the superposition of the linear equilibrium trend and a single cosine, which may mathematically represent local changes in bed elevation induced by large-scale bedforms: where A and L indicate initial wave height and initial wave half-length, respectively. From Equation (26) with no source term, one obtains: where N = 0 (that is, at x = x 0 = L/2), section-averaged bed elevation coincides with the elevation obtainable from the only reach-scale trend: η(x 0 , t) = a − i f eq x 0 . For x < x 0 , N-profile is convex: ∂ 2 N/∂x 2 < 0 and the small-scale three-dimensional undulations represented by η migrate upstream (v(x, t) < 0) with an exponentially decaying velocity; conversely, for x > x 0 , N-profile is concave: ∂ 2 N/∂x 2 > 0 and the small-scale undulations migrate downstream (v(x, t) > 0) at the same gradually decreasing rate. Note that the convex N-profile is associated with N > 0 and η(x, t) > a − i f eq x; the vice versa holds for the concave counterpart. Considering that the critical state of flow implies i f eq ≡ i c , where i c indicates reach-scale critical bed slope, the quantity η c = a − i f eq x may be assimilated to the reach-scale critical bed elevation. In this case, hinging on basic hydraulics (e.g., [51]), the condition η(x, t) > η c (i.e., the condition of local section-averaged bed elevation larger than the critical one) can be associated with a locally slow (sub-critical) sediment wave. The vice versa can be said in the case of concave N-profile, which would turn out to be associated with a locally fast (supercritical) sediment wave ( Figure 1). Now, again from basic hydraulics, we know that small perturbations (here represented by η ) typically go upstream in the presence of sub-critical free-surface flows; conversely, they go downstream in the presence of supercritical flows. Therefore, that result makes our model perfectly compatible with a sediment-stream version of the classic gradually-varied flows theory. Specifically, we know that the approximate integration of de Saint-Venant equations for rectangular channel section [57]: for i f ∼ = J leads to: and U ± gH ∂H ∂x Note that Equation (38) represents free-surface waves that propagate with celerity: Therefore, for slow flows (U/ gH < 1) and dU/dH < 0, c < 0 and the perturbations go upstream; on the other hand, for fast flows (U/ gH > 1) and dU/dH < 0, c > 0, and the perturbations go downstream.

Results
To solve system (26)- (27) showing that it can explain the tendency of real rivers to maintain braiding under particular conditions and follow and monitor their evolution from an equilibrium configuration to the other, we resorted to a mixed analytical-numerical method. Specifically, we solved Equation (26) analytically for S(x,t) = 0 and an initial condition represented by an oscillatory pseudo-periodic function that is qualitatively consistent with our laboratory experimental observations [3]: where ∆ indicates the amplitude and Λ the characteristic length-scale, obtaining: Then, we solved Equation (27) by particle-tracking (e.g., [58][59][60][61]) for s (x, y, t) = 0 and: where ξ has the dimensions of a length, and h(x) and m(y) are two numerical coefficients defined as: The reference macro-data set belongs to a river Basento reach (Basilicata region, southern Italy, Figures 2 and 3) that can be considered as truly braided according to several classification criteria [62][63][64][65][66]   Additionally, we assumed ∆ = 1 m, ξ = 2d = 0.09 m, and β = d/2 = 0.0225 m. Therefore, the initial bed elevation deviation turned out to be represented by a 9 cm-high double array of sediment particles instantaneously placed about the channel axis between x = 0 and x = B = 150 m (total number of particles: 6667), with a consistently almost zero average deviation: Provided that the linear structure of (26) and (27) enables the superposition principle, the aim of the present numerical experiment was to investigate how an elementary streamwise perturbation evolves toward different types of braided or unbraided configurations as a function of the crucial ratio Λ/B. The application to real initial conditions will be the object of future research, when evolving post reach-equilibrium DEMs (Digital Elevation Models) will be available from the literature or our ongoing supplementary laboratory experiments.
For particle-tracking purposes, the barycenter of each sediment particle was virtually tracked based on the following discretized solution of the trajectory equation: where X and Y, respectively, indicate the stream-wise and transverse component of the particle position vector, G (0,1) indicates a zero-mean and unit-variance Gaussian distribution, dt is the elementary time step and: v( The elementary time step was chosen in such a way as to limit the length of the elementary purely diffusive displacement to a small percentage of channel total width: with: When a particle overcame a side boundary (representing the impervious bank), it was back-reflected by |dY| = √ 2Ddt 2 G(0, 1), where dt 2 = dt − dt 1 and dt 1 was the time needed to reach the boundary itself. Note that, during dt 2 , the particle was not allowed to move along the longitudinal direction. It should be emphasized that, although our particletracking is 2-D in the horizontal projection of the channel plan view, the effect of lift and drag is inherent in the velocity and diffusion coefficient of the advection-diffusion equation that it solves. There are in the literature a few examples of 3D particle tracking applied to the study of saltating particles diffusion in turbulent water flow in more general conditions (e.g., [67]). The peculiarity of our study consists in the specific and more restricted field of application (braided gravel beds) and in the analytical nature of the model, allowed by the combination of the high grain Reynolds numbers quadratic relationship between dimensionless unit bedload rate and stream power and the typical critical state of flow, which determines the way the stream power is evaluated. Figure 4 shows the velocity distribution sampled for Λ = B/8 by the double sediment array at t = t 0 and t = τ d , with: representing a characteristic diffusive time (i.e., the time needed by pure diffusion to cover a distance proportional to Λ). Figure 5 shows the solution of Equation (27) in terms of particles' barycenter position at t = τ d . Finally, Figure 6 shows the corresponding Xhistogram. As one can see, the high-frequency oscillations of velocity induce a marked fragmentation of the initial sediment strip and a relatively reduced diffusion. The histogram displaying X-frequency (f ) distribution exhibits an almost two-peaked shape with a central minimum. Figures 7-9 show the same as Figures 4-6 but referred to Λ = B/4. In this case, the initial velocity was less oscillating and overall characterized by a smaller magnitude. The number of sediment accumulation sites was reduced while their dimensions were increased due to the gradual transition from advective to diffusive regime. The X-histogram revealed a multi-peaked structure with higher concentrations of sediment in a lower number of longitudinal bands. Figures 10-12 refer to Λ = B/2. The initial velocity was now rather low and only weakly oscillating. The corresponding particles' distribution after the actual τ d highlights the formation of three well-defined sediment accumulation sites in about 190 m, associated with three histogram peak zones. Note that the case-study in Figure 3, which represents a clearly braided branch of river Basento, and was characterized by Λ/B ∼ = 1/2, include three large bars in about 220 m. The case Λ = B (Figures 13-15) was the first characterized by a monotonic initially sampled velocity. From Figure 14, one can appreciate the dominant effect of diffusion (many sediment particles had already reached the boundaries at t = τ d ) with considerable back-displacements (X < 0) and the slow transformation of the bed configuration into a two-dimensional sand-waves pattern. This tendency finds its maximum expression for Λ = 2B (Figures 16-18), with channel width that is now almost uniformly occupied by sediment particles and the X-histogram that is rather close to a (skewed) Gaussian distribution, as in Taylor-like diffusion. Since the numerical experiments were conducted with the same B, Q, i feq , ρ s , p, and d and, therefore, with the same diffusion coefficient D, we can conclude that the tendency to maintain braiding and flow bifurcation is associated with equilibrium average bed profiles N and, therefore, equilibrium average stream power γQ i f eq − ∂N/∂x characterized by maximum period Λ that does not exceed transverse channel dimension.               Finally, in Table 2, we summarize the main numerical parameters and results: the ratio Λ/B, the diffusion time τ d , the Péclet number Pe = Λv max /D-which is representative of the relative importance of advection and diffusion-and the braiding index BI, which was defined as the ratio of channel length occupied by sediment particles at t = τ d to velocity characteristic length Λ. As one can see, the smaller the ratio Λ/B, the smaller the time needed for diffusion over Λ (τ d ), the larger Pe and the larger BI. For the sake of completeness, and to analyze the free evolution of the equilibrium section-averaged bed profile,  show N(x) at t = t 0 and t = τ d for Λ = B/8, Λ = B/4, Λ = B/2, Λ = B and Λ = 2B, respectively. As one can see, due to the specific structure of (26), the pseudo-periodic behavior of the initial distribution tends to approach a perfectly periodic (other than considerably attenuated and characterized by lower frequency) oscillation, which is in turn affecting the local deviation evolution shown in the previous figures.

Discussion and Conclusions
In the present study, we propose an analytical-numerical methodology for the reconstruction of the equilibrium dynamics of gravel braided beds. The mathematical model is based on the combination of sediment continuity with or without forcing term and the practically quadratic relationship between dimensionless unit bedload rate and dimensionless unit stream power, which was derived by the authors in a previous work by regression, including ad hoc laboratory experimental measurements performed at high grain Reynolds numbers and real river data available in the literature. The crucial assumption enabling a semi-analytical solution in terms of bed elevation function was represented by the recognition of the critical state of flow, widely documented in the literature, which makes local section-averaged energy loss coinciding with local section-averaged bed slope. The first important result of the study is represented by the specific structure of the bed elevation deviation governing equation, coupled with a variable-coefficient diffusion equation holding for the section-average component. Note that the diffusive nature of sediment transport had already been highlighted by previous experimental investigations. A partial simplification of the two coupled differential equations was achieved in terms of diffusion coefficient, which was found to be proportional to the sum of reach-scale bed slope, section-averaged bed slope, and deviatory bed slope. As a matter of fact, due to the pseudo-periodical oscillation of section-averaged bed elevation and corresponding deviation, which was experimentally verified, the diffusion coefficients practically include terms that determine a mutual compensation of mass transfer from higher to lower concentrations of sediment and vice versa. For that reason, the dependence of these diffusion coefficients on the oscillatory part of the bed elevation function was neglected, and the only dependence on the constant reach-scale bed slope was maintained. A second important result, which basically proves the consistency of the present mathematical formulation, is represented by the direction of the small-scale bed forms migration according to the local section-averaged bed curvature. Specifically, where the section-averaged bed is convex, sub-critical bed waves would move upstream; conversely, where the section-averaged bed is concave, super-critical bed waves would migrate downstream. This conclusion is fully compatible with classic hydraulics and gradually-varied flows theory. The system of coupled equations that govern section-averaged bed elevation and corresponding deviation was analytically-numerically solved for a set of macro-data belonging to a real truly braided river reach of southern Italy, assuming an initially pseudo-periodical sectionaveraged bed. The advection-diffusion equation governing the bed elevation deviation was solved by particle-tracking, based on the constant diffusion coefficient and the velocity values locally sampled by each sediment particle, with an initial condition represented by a central double array of stones. The solution was plotted in terms of particles' barycenter position and corresponding histogram after a characteristic diffusive time, for five different values of the ratio of (initial) longitudinal macro-bedforms maximum period to channel maximum width. The results indicate that the smaller the ratio, the more accentuated the tendency of the river to maintain the braided bed structure, with the formation of smaller and highly concentrated sediment heaps whose number is proportional to the frequency of the oscillatory section-averaged bed profile. Conversely, when the ratio is large (specifically, larger than 1), bed dynamics tends to approach a Taylor-like dispersive process (e.g., [68]) with asymmetric tails (e.g., [69]). This result means that the initial disturbance being the same, the longitudinally smoother the equilibrium braided bed (longer macro-bedforms or bars), the more marked the diffusive nature of transport, with the displacement of particles in all directions that are only bounded by river valley side boundaries. Conversely, the longitudinally more heterogeneous the equilibrium section-averaged bed (shorter macro-bedforms or bars), the more marked the advective nature of transport, with mainly longitudinal particle displacement and sediment mass fragmentation. The physical explanation for this behavior is that lower-frequency section-averaged bed waves are associated with lower-amplitude advective particle displacements and, therefore, relatively more intense diffusion. Finally, it should be noted that the system of coupled differential equations derived in the present study may easily be utilized (eventually in its complete version (23)-(24)) for a finite-differences/elements solution, starting from the real DEM of the river bed at t = 0. As a matter of fact, it is widely recognized that human-induced factors, such as embankments, roads, bridges, and dams construction, can affect sediment transport and bed morphology of the regulated river [70][71][72], and can critically impact on location and extent of gravel bars [73], which are important natural structures for fluvial or riparian habitats and, therefore, for the protection and conservation of species of fishes, invertebrates, plants, and birds [74][75][76]. Several fish species live in braided rivers, and their connected wetlands and gravel bed grounds constitute suitable spawning areas. Additionally, in many cases, gravel bars represent an important substrate for the establishment and development of ground flora and woody vegetation and guarantee higher plant diversity [77]. Sustainable management of braided rivers should, therefore, ensure their ecological potential and biodiversity by preserving braiding structure over time [1]. The analytical-numerical solution proposed by the present study, which models gravel braided bed evolution from an equilibrium configuration to the other for imposed flow and bedload rate, may, therefore, be applied for the analysis of the impact of fluvial engineering works on river habitat and its biotic components, as well as for the planning and the monitoring of sustainable river restoration.