A One-Way Coupled Hydrodynamic Advection-Diffusion Model to Simulate Congested Large Wood Transport

: An advection-diffusion model is proposed to simulate large wood transport during high ﬂows. The mathematical model is derived from the wood mass balance, taking into consideration both the wood mass concentration and the log orientation, which affects log transport and, most importantly, wood accumulation. Focusing on wood mass transport, the advection-diffusion equation is implemented in a hydrodynamic model to provide a one-way coupled solution of the ﬂow and of the ﬂoating wood mass. The model is tested on a large series of ﬂume experiments, involving at least 30 logs and different control parameters (ﬂow Froude number, log length, diameter, release point). The validation through the experimental data shows that the proposed model can predict the correct displacement of the most probable position of the logs and to simulate with a sufﬁcient accuracy the planar diffusion of the wooden mass. Transversal wood distribution is more accurate than the streamwise one, indicating that a higher control on the longitudinal diffusion needs to be implemented.


Introduction
Rivers are important components of urban settings. In the past, human relied on rivers for protection and livelihood while presently their role is mainly related to social and ecological benefits, through several revitalization projects all around Europe [1]. Hydraulic risk is a counterpart of the proximity of urban and river environments. Channel narrowing, due to bridges or culverts, represents a critical section in case of flood [2], while the increase of impervious areas [3] and hydraulic implications of climate change, either increased water discharge or increased number of extreme floods [4][5][6], may affect flood hazard, requiring prompt adaptation or redaction of management plans. Governments and supernational groups urge to determine areas at risk and countermeasures to protect lives and material assets, to reduce losses and associated costs.
Flood hazard maps are a fundamental tool for this purpose, since they help in identifying which areas are at risk and to what extent planned measures are effective. They also represent a powerful tool for the communication of risk [7,8], to raise public awareness about endangered areas and proper behavior in case of flooding.
The adoption of two-dimensional modeling and the use of high-resolution Digital Terrain Model enable an accurate computation of water depth and velocity both in the channel and in floodplain areas [9]. To further improve the accuracy of flood risk estimation, additional issues that occur during floods should be included in hydraulic risk modeling, such as the transport of sediments, pollutants, or wood.

Advection-Diffusion Model for Large Wood
During high flows, the congested transport of large wood involves the motion of many logs and wooden debris, behaving like a continuum wood carpet that floats on the water surface. The total mass in an infinitesimal planar area dxdy is the product of wood density ρ w and the probability density function p(x,θ,t), where the state-space representation includes the planar coordinates x = (x,y) and the planar orientation θ, besides time t. The planar orientation of logs is included since it affects the interaction between logs and with the riverbanks and, consequently, the actual wood mass in the infinitesimal area. More logs fit the same area, if they are aligned with the streamlines, than in cross-flow configuration.
The wood mass continuity equation results in the form of the advection-diffusion equation: where p(x,θ,t) is the probability density function of wood in the space (x,θ), v is the planar transport velocity, ω is the angular velocity, K x and K y are the diffusion coefficients in the space (x,y) for anisotropic diffusion, K θ is the angular diffusion coefficient and S w represents a source/sink term for the wood probability density function, like a source of mass in the domain that does not affect the flow field, since, as previously stated, wood is considered to be a passive transported substance, one-way coupled with the water flow.
Since two movements occur during wood transport, translation and rotation, the probability density function p(x,θ,t) can be considered to be the product of two distinct contributions: where c w (x,t) refers to the probability that a wooden element occupies the position x, independently on its orientation, while δ(x,θ,t) gives the probability that a wooden element, positioned in x, has an orientation ranging in (θ,θ + dθ). Substituting Equation (2) into Equation (1) and neglecting the source term S w , we obtain: By integrating Equation (3) in the interval (0,2π), being 2π 0 δ(x, θ, t)dθ = 1, due to the flow continuity equation and to the fact that c w (x,t) and v are independent of θ, it results: which is the advection-diffusion equation of the positional probability density function. Finally, multiplying Equation (4) by δ, subtracting it to Equation (3), and dividing the results by c w (x,t), the advection-diffusion equation of the angular probability density function reads: At this stage, the main interest is in the transport of wooden mass, since we aim at evaluating the capabilities of the model to replicate wood transport and not the interactions with inline structures. For this reason, the solution of Equation (5) is not considered in this paper, although its implementation is fundamental for future steps and for the simulation of wood accumulation at bridge piers.

Outlines of the Coupled System
To provide an accurate solution of the flow and LW dynamics, the advection-diffusion model for mass probability density function (Equation (4)) is coupled with the solution of the SWE. The coupling is performed following the mathematical model proposed in [33]. Here, for briefness, only the final form of the system of equation that needs to be solved is reported: ∂U ∂t where h is the water depth, q x = uh and q y = vh are the unit discharges along the x and y coordinates, u and v are the components of the velocity vector v, along the x and y coordinates, respectively, g is the gravity acceleration, S 0 is the source term due to bed slope, S f is the source term due to friction (expressed with the Manning coefficient), K is the diffusion matrix (diagonal matrix) and ϕ is the volume averaged wood mass. The volume averaged wood mass is obtained by multiplying the spatial probability density function c w (x,t) by the wood density ρ w and dividing it by the water volume of the computational cell: where dxdy represents an ideal planar area. The numerical model ORSA2D is thus modified to solve the conservative form of Equation (6). ORSA2D is a hydrodynamic model that follows a finite-volume approach, implementing a Roe-Riemann scheme to solve the SWE [34]. It can model both steady and unsteady flow, with proper boundary conditions that can be constant or time dependent. The numerical approach for the coupled solution is essentially the same adopted for the solution of the SWE alone [34], the main difference being the additional unknown included in the new system of Equations (6), which results in an additional eigenvalue and an additional eigenvector, plus one element for the others, of the approximated flux Jacobian.
The details of the numerical model are the same used by [33] in the coupling between SWE and solute flow and are fully described there. Here, we only report the resulting approximated eigenvalues and eigenvectors for the coupled SWE-LW case: In Equations (9) and (10), the eigenvalue λ 4 refers to the wood transport equation, as well as the fourth term of the eigenvectors e 1 , e 2 , e 3 and the eigenvector e 4 . The symbols with tilde used in Equations (9) and (10) are the averaged values of the velocity components (u,v), celerity c and wood averaged concentration ϕ: where the R and L subscripts refer to the considered side of an edge of the computational cell (respectively, right, or outer, and left, or inner, side of the edge, Figure 1).
where the R and L subscripts refer to the considered side of an edge of the cell (respectively, right, or outer, and left, or inner, side of the edge, Figure The boundary conditions for the system of equation Equation (6), cons are the standard ones usually included in hydrodynamic models (e.g., wat charge hydrograph, rating curve, Froude number hydrograph [34]).
Regarding wood, a wood mass hydrograph can be set as upstream bo tion while a simplified condition is proposed for the downstream bound concentration at the right side of a boundary edge (φR) is set equal to the t centration of the computational cell, so that the average value computed (Equation (11)) depends only on the difference between the left and righ Since, for boundary edges, the right water depth hR is calculated from th boundary conditions, the wood mass follows the variations of the water related to the selected condition.
The diffusion source term in Equation (6) is treated as in [35], integr computational cell, applying Gauss's theorem and then discretizing the co with a sum over the cell edges: where n is the normal vector of the edge (pointing outside the computatio the number of edges, s is the length of each edge, KR and KL are the matrix diffusion, h * R and h * L are the water depth (predictor values) at each edge s φL) is the difference between the volume averaged wood mass (evaluated w method) and dRL is the distance between the center of the cells right and le to the considered edge. R and L subscripts are the same as above.
In addition, to reduce the numerical diffusion observed for the wood and improve the solution accuracy of the diffusion term also with a first numerical scheme, a numerical correction was derived from [35]. It involv tation of a corrective diffusion matrix Kcorr, which is only dependent on the m and on the fluid velocity that is subtracted to the diffusion term if it is positi than it. The correction is calculated for each time step Δt: and it is implemented as follows: The boundary conditions for the system of equation Equation (6), considering water, are the standard ones usually included in hydrodynamic models (e.g., water stage or discharge hydrograph, rating curve, Froude number hydrograph [34]).
Regarding wood, a wood mass hydrograph can be set as upstream boundary condition while a simplified condition is proposed for the downstream boundary. The wood concentration at the right side of a boundary edge (ϕ R ) is set equal to the total wood concentration of the computational cell, so that the average value computed for calculation (Equation (11)) depends only on the difference between the left and right water depth. Since, for boundary edges, the right water depth h R is calculated from the downstream boundary conditions, the wood mass follows the variations of the water depth that are related to the selected condition.
The diffusion source term in Equation (6) is treated as in [35], integrating it on the computational cell, applying Gauss's theorem and then discretizing the contour integral with a sum over the cell edges: where n is the normal vector of the edge (pointing outside the computational cell), NE is the number of edges, s is the length of each edge, K R and K L are the matrix for anisotropic diffusion, h R and h L are the water depth (predictor values) at each edge side, δϕ = (ϕ R − ϕ L ) is the difference between the volume averaged wood mass (evaluated with an explicit method) and d RL is the distance between the center of the cells right and left, with respect to the considered edge. R and L subscripts are the same as above.
In addition, to reduce the numerical diffusion observed for the wood mass transport and improve the solution accuracy of the diffusion term also with a first order accurate numerical scheme, a numerical correction was derived from [35]. It involves the computation of a corrective diffusion matrix K corr , which is only dependent on the mesh geometry and on the fluid velocity that is subtracted to the diffusion term if it is positive and smaller than it. The correction is calculated for each time step ∆t: and it is implemented as follows: Finally, the time step is calculated considering both the Peclet and the Courant-Friedrichs-Lewy (CFL) numbers, limiting the value to further reduce diffusivity problems [35].
The diffusion coefficients of the matrices K R and K L are derived from the results presented in [30], in which streamwise and transversal coefficients (K s , K t ) were provided as a function of the particle Froude number Fr p and of the transversal relative release distance t Rr of the mass concentration (Equation (16)). These values are combined according to the local velocity orientation to obtain the elements of the diffusion matrices (Equation (17)). where , L w is the length of the wood pieces, t R is the release distance from the right bank, u and v are the flow velocity component, K xx and K yy are the diagonal element of the diffusion matrix and K xy the off-diagonal elements.

Experiments Description
The numerical model is tested on a series of flume experiments performed at the University of Trento, described in [30]. Here the main aspects relevant to the comparison with the numerical results are reported.
Two flume configurations were employed: a S-shaped flume and a S-shaped flume with a gradual Venturi narrowing that halves the flume section. In the first case, the bottom slope was also varied. The flumes are 2 m wide, with a reduction to 1 m in the narrowing, and about 22 m long, with rectangular section, vertical side walls (1 m high) and a fixed gravel bed (D 50 = 8 mm).
The experiments were performed in steady flow conditions, releasing logs one by one and filming their trajectories with 3 GoPRO Hero 5 cams, located 2.5 m above the flume (30 fps, pixel resolution 1920 × 1080, linear Field Of Vision), able to frame altogether nearly the entire length of the flume. Figure 2 shows the flume planar views and a frame from the central camera. The flow velocity was measured with an Acoustic Doppler Velocimeter for each hydraulic condition, in the points shown in circles in Figure 2.  Through data analysis [30] each log trajectory was estimated, acquiring time, planar positioning, and orientation. Cylindrical logs, smooth and varnished with green impermeable painting, were employed. Lengths and diameters were selected based on a theoretical physical scale (λ h = 60 and λ h = 40, horizontal and vertical, respectively). The cylinders average density is 790 kg m −3 , such as that of green wood broadleaved trees [36]. Table 1 resumes the characteristics of the experiments used as a reference for the numerical model. The main parameters relevant to the experimental analysis and numerical simulation are the flow Froude number Fr = v mean √ gh , which remains constant only for the S-shaped flume, the transversal position of release (calculated from the right sidewall) t R , the log length L w and diameter D w . The flow velocity was measured with an Acoustic Doppler Velocimeter for each hydraulic condition, in the points shown in circles in Figure 2.

Modeling the Transport Velocity
In the considered experiments, and in most cases in the natural environment, wood density is lower than the water density [30], so wood logs float on the water surface and are transported by the superficial flow. Under such conditions, the vertical averaged mean velocity, computed by the SWE model, is not representative of the actual transport velocity v of the advection-diffusion model (Equation (4)). This velocity is more likely to be the maximum flow velocity, computed as a function of the mean value assuming, for example, a (1/7)th power velocity distribution: In addition, when a log is entrained in the flow, its initial velocity is lower and tends to adapt to the surface flow velocity. The same behavior was observed in the experiments described in Section 2.2, which presented a velocity near to 0 m s −1 at the beginning due to the permanence in the release device. In [30], an asymptotical formulation is proposed based on the analysis of a single test case (Fr = 0.4, L W = 0.20 m, D W = 0.02 m, t R = 0.80 m). To generalize that formulation, a non-dimensional coefficient c vel is obtained: that can be multiplied by the mean flow velocity of each computation cell containing wood mass concentration to estimate the actual transport velocity (v = c vel v mean ). The value c at the exponent denominator is linearly correlated with the Froude number Fr R , evaluated at the release point, and its expression is: Logs velocity estimation in the S-shaped flume with narrowing showed that just downstream of the narrowing, the logs tend to maintain a velocity higher than the theoretical one (computed with the asymptotical formulation). This probably occurs because of the log inertia that delays the adaptation of the log velocity in case of sudden flow deceleration, especially in presence of section variation, energy loss and turbulent flow behavior. An additional coefficient is thus included in the proposed model to increase wood velocity in case of local and abrupt decelerations, i.e., when the maximum wood concentration undergoes a velocity reduction greater than a threshold value based on the experiments. This coefficient is related to the release transversal position t R : if the logs are nearer to the flume sidewalls, they are subject to higher gradients that tend to diminish their velocity, while if they flow nearer to the flume axis, they maintain a higher velocity. The deceleration coefficient k d is thus defined as follows: Overall, the transport velocity computed at each time step, for each cell i containing a wood concentration, is: where Fr R is the Froude number at the release point, v f is the flow velocity computed by ORSA2D, i c max is the cell where the maximum concentration is found at the considered time instant and ∆v is a threshold value (∆v = −0.16 m s −1 based on the experimental tests).
The transport velocity resulting from Equation (22) is shown in Figure 3, compared with the experimental log velocity, with the mean flow velocity computed by ORSA2D and with the maximum theoretical flow velocity obtained by applying Equation (18).

Details of the Numerical Tests
The experiments described in Section 2.2 are simulated with ORSA2D model modified to simulate wood mass advection-diffusion. Both the S-shaped flume and the Sshaped flume with narrowing are discretized with triangular cells of about 0.05 m side, for a total number of 39,856 elements for the first, and 38,256 for latter (which presents fewer elements since the narrowing is treated as a solid physical boundary, see Figure 4).

Details of the Numerical Tests
The experiments described in Section 2.2 are simulated with ORSA2D model modified to simulate wood mass advection-diffusion. Both the S-shaped flume and the S-shaped flume with narrowing are discretized with triangular cells of about 0.05 m side, for a total number of 39,856 elements for the first, and 38,256 for latter (which presents fewer elements since the narrowing is treated as a solid physical boundary, see Figure 4). The inlet discharge, set as upstream boundary condition, and the bottom slope are shown in Table 1. The roughness is represented with a constant Manning coefficient, n = 0.021 s m −1/3 , while the imposed downstream boundary condition is the downstream Froude number, set equal to 0.09, 0.25 and 1 for the three tests in the S-shaped flume, and 0.36 for the flume with narrowing. While the latter was estimated from experimental measures, the first three were calibrated to reduce the differences between the experimental and the numerical flow velocities.
The wood mass is assumed to be released in a single cell, corresponding to the mean initial position of the logs for each test, with a triangular wave function. Equation (23) shows the computation of the input concentration rate: where the numerator represents the total mass for the considered experiment (NL being the number of repetitions, or of total logs released, as in Table 1) while the denominator includes the input cell area, Acell, water depth hcell and the estimated release duration Tref. The release duration is set equal to 1.2 s and was evaluated based on the actual release time observed during the experiments performed for the experimental campaign of [32].

Hydraulic Simulation
Numerical results are compared with the available measures (velocity measures in the points shown in Figure 2) to provide an evaluation of the simulation accuracy. Water depth was not measured for the S-shaped flume, so the local Fr was not computed.
Regarding the S-shaped flume ( Figure 5), a good agreement is observed for the measured and simulated velocities. Better correlation is observed for the v component ( Figure  5b, average correlation coefficient R 2 = 0.71) than for u component (Figure 5a, average correlation coefficient R 2 = 0.20), which is less variable in the simulation than in the experiments (especially for the lower Froude numbers). The inlet discharge, set as upstream boundary condition, and the bottom slope are shown in Table 1. The roughness is represented with a constant Manning coefficient, n = 0.021 s m −1/3 , while the imposed downstream boundary condition is the downstream Froude number, set equal to 0.09, 0.25 and 1 for the three tests in the S-shaped flume, and 0.36 for the flume with narrowing. While the latter was estimated from experimental measures, the first three were calibrated to reduce the differences between the experimental and the numerical flow velocities.
The wood mass is assumed to be released in a single cell, corresponding to the mean initial position of the logs for each test, with a triangular wave function. Equation (23) shows the computation of the input concentration rate: where the numerator represents the total mass for the considered experiment (N L being the number of repetitions, or of total logs released, as in Table 1) while the denominator includes the input cell area, A cell , water depth h cell and the estimated release duration T ref .
The release duration is set equal to 1.2 s and was evaluated based on the actual release time observed during the experiments performed for the experimental campaign of [32].

Hydraulic Simulation
Numerical results are compared with the available measures (velocity measures in the points shown in Figure 2) to provide an evaluation of the simulation accuracy. Water depth was not measured for the S-shaped flume, so the local Fr was not computed.
Regarding the S-shaped flume ( Figure 5), a good agreement is observed for the measured and simulated velocities. Better correlation is observed for the v component (Figure 5b, average correlation coefficient R 2 = 0.71) than for u component (Figure 5a, average correlation coefficient R 2 = 0.20), which is less variable in the simulation than in the experiments (especially for the lower Froude numbers).

Wood Transport Simulation
The wood transport simulations are carried out from steady state hydraulic conditions. The wood concentration is released at 50 s, which is then considered to be the initial time for wood transport. As shown in Table 1, several tests were modeled that differ for the main parameters affecting diffusion: flow Froude number, log length L w , log diameter D w and release distance t R from the right side of the channel.
For each hydraulic condition examined in the previous section, the results of selected tests are presented, at regular time intervals and at the last available instant, while in the discussion paragraph a synoptic comparison will be shown. Please note that in the experiments each log was separately released, following its own trajectory, and not interacting with the others. Grouped representations have the sole purpose to help the distribution visualization and results discussion.

Wood Transport Simulation
The wood transport simulations are carried out from steady state hydraulic conditions. The wood concentration is released at 50 s, which is then considered to be the initial time for wood transport. As shown in Table 1, several tests were modeled that differ for the main parameters affecting diffusion: flow Froude number, log length Lw, log diameter Dw and release distance tR from the right side of the channel.
For each hydraulic condition examined in the previous section, the results of selected tests are presented, at regular time intervals and at the last available instant, while in the discussion paragraph a synoptic comparison will be shown. Please note that in the experiments each log was separately released, following its own trajectory, and not interacting with the others. Grouped representations have the sole purpose to help the distribution visualization and results discussion.

S-Shaped Channel with Venturi Narrowing
Graphs in Figure 8 show the comparison of the experimental log position and the computed wood concentration for three tests with different release distance from the right bank, in the flume with the narrowing. The time sequence has an interval of 5 s and the minimum value of the volume averaged concentration that corresponds to the largest contour line, is 5 kg m −3 . The final time, which will be considered to be a reference time for the discussion, varies slightly for each test (71 s, 72 s and 74 s).

S-Shaped Channel with Venturi Narrowing
Graphs in Figure 8 show the comparison of the experimental log position and the computed wood concentration for three tests with different release distance from the right bank, in the flume with the narrowing. The time sequence has an interval of 5 s and the minimum value of the volume averaged concentration that corresponds to the largest contour line, is 5 kg m −3 . The final time, which will be considered to be a reference time for the discussion, varies slightly for each test (71 s, 72 s and 74 s).  Figures 7 and 8 show how the wood is transported in the experimental tests. It must be highlighted that despite the experiments were performed releasing the logs one by one [30], in these figures they are shown altogether as black segments, to better visualize their distribution. In general, logs are grouped at the initial time, since they are released in the  Figures 7 and 8 show how the wood is transported in the experimental tests. It must be highlighted that despite the experiments were performed releasing the logs one by one [30], in these figures they are shown altogether as black segments, to better visualize their distribution. In general, logs are grouped at the initial time, since they are released in the same position, and then spread longitudinally and transversally because of flow advection and surface turbulent diffusion.

Discussion
Similarly, the numerical results are concentrated at the initial time and tend to expand while flowing downstream. In all the simulations, the planar enlargement of the area of non-zero concentration (red contours in Figures 7 and 8) appears to be wider than the scattering of the experimental logs, although the outer contour line represents a very low volume averaged concentration (5 kg m −3 ) that for the volume of an average cell, corresponds to a wood mass lower than 1 g.
For the S-shaped flume (Figure 7), the peak of the numerical concentration often moves downstream faster than the experimental logs, as it can be viewed by the center of the contour lines being located downstream of the center of mass of the groups of logs in the experiments. The opposite happens in case of constriction (Figure 8), where at each time step, a better agreement between the area occupied by the logs and by the numerical contour lines is found.
To provide a synoptic evaluation of the model, a comparison between the experimental concentration and the numerical tests at the final time instant is considered. The time is selected as the maximum time at which all the logs released in the flume (N L ) are visible in the recording (e.g., the final time shown in the Figures 7 and 8, but it varies for each test).
To analyze the accuracy of the simulation, the location of the mean position of the experimental logs is compared with the location of the cell with the maximum concentration (i.e., the most probable position) in the simulation. Figures 9a and 10a show that the mean planar position for each test is, on average, in the range 14.    The total mass in the target area is shown in Figure 11. For both the flumes, the simulated mass corresponds to the experimental one within a confidence interval of 40%. In all cases, the total simulated mass is equal or lower than the experimental one. Please note that the total mass at the end of each simulation is equal to the released mass, so the reduction observed in Figure 11 is related to the dimension of the target area, which is too small to comprise the total concentration. The extension of the target area, in fact, is set upon the experimental planar distribution, which is clearly reduced with respect to the  (Figure 10b). Considering both flumes, the global mean absolute percentage error is 4.5% and 16.5% for the streamwise and transversal coordinates, respectively.
Up to now the focus was on the log mean position, or most probable location of the maximum concentration. An additional validation can be performed by checking if the model can distribute the mass as observed in the experiments.
For this purpose, a "target area" is identified, as a series of concentric rings, which are centered on the log mean position for the experiments and on the coordinates of the cell with the maximum concentration for the simulations. The rings are 0.05 m wide and the maximum diameter is set to include the maximum span of the log positioning (which may differ for each test).
For the experiments, the concentration is computed by dividing the total mass in each ring (number of logs inside the ring times log volume times wood density, M = N in ring L W D W 2 4 ρ W ) by the ring planar area. The value obtained is the planar concentration of wood mass, named "experimental concentration" hereafter. For the numerical simulations, the total concentration in each ring is the sum of the concentration of each domain cell that is included in that concentric ring.
The total mass in the target area is shown in Figure 11. For both the flumes, the simulated mass corresponds to the experimental one within a confidence interval of 40%. In all cases, the total simulated mass is equal or lower than the experimental one. Please note that the total mass at the end of each simulation is equal to the released mass, so the reduction observed in Figure 11 is related to the dimension of the target area, which is too small to comprise the total concentration. The extension of the target area, in fact, is set upon the experimental planar distribution, which is clearly reduced with respect to the numerical one (e.g., see Figures 7 and 8). Independently of the flume configuration, the lowest accuracy is obtained for the experiments with the largest mass involved (i.e., wooden samples with L W = 0.4 m and D W = 0.03 m). On average, the numerical wooden mass included in the target area is about 85% of the experimental one.   Figures 12 and 13 show the difference of the numerical concentration m perimental one, in each ring. In most cases the largest differences are in the with the numerical values being lower than the experimental ones for both fl urations. Once again, the simulations diffuse the mass concentration more th with lower values in the central rings and a concentration distributed on an than the target one.
The position of the centers of mass of the logs in the experiments, sho dots in all the graphs of Figures 12 and 13, confirms that the logs are not evenly in the circular areas but present, in general, a longitudinal dispersion high transversal one, which is more or less evident depending on the considered contrary, the numerical results are more gradually distributed. As well as mass estimation, the largest differences are observed for the tests involvin logs, so the greater mass.
A quantitative evaluation of the difference of the distributed planar con computed with the Mean Absolute Error (MAE = |exp.conc. − num.conc.|, valua concentric ring), which has a global average value of 0.22 kg m −2 .
Finally, the wood mass distribution along the two main diameters of th (in streamwise and transversal direction) is shown in Figure 14. All the lines on the position of the maximum mass, in the streamwise or transversal direc slot amplitude to calculate the mass in each spatial interval (either sum of th or sum of the mass in the cells) is of 0.05 m. The maximum width of the g Figure 11. Comparison of the experimental and numerical mass included in the target areas for the two flume configurations. Figures 12 and 13 show the difference of the numerical concentration minus the experimental one, in each ring. In most cases the largest differences are in the inner rings, with the numerical values being lower than the experimental ones for both flume configurations. Once again, the simulations diffuse the mass concentration more than expected, with lower values in the central rings and a concentration distributed on an area wider than the target one.
The position of the centers of mass of the logs in the experiments, shown as black dots in all the graphs of Figures 12 and 13, confirms that the logs are not evenly distributed in the circular areas but present, in general, a longitudinal dispersion higher than the transversal one, which is more or less evident depending on the considered test. On the contrary, the numerical results are more gradually distributed. As well as for the total mass estimation, the largest differences are observed for the tests involving the longer logs, so the greater mass.
A quantitative evaluation of the difference of the distributed planar concentration is computed with the Mean Absolute Error (MAE = |exp.conc. − num.conc.|, valuated for each concentric ring), which has a global average value of 0.22 kg m −2 .
Finally, the wood mass distribution along the two main diameters of the target area (in streamwise and transversal direction) is shown in Figure 14. All the lines are centered on the position of the maximum mass, in the streamwise or transversal direction, and the slot amplitude to calculate the mass in each spatial interval (either sum of the logs mass or sum of the mass in the cells) is of 0.05 m. The maximum width of the graphs corresponds to the maximum radius of the target area. This suggests that the model provides an acceptable simulation of the transversal diffusion, while the longitudinal one is overestimated.  Table 1 for the order of the experiments.  Table 1 for the order of the experiments. Figure 14 confirms that independently of the flume shape, the experimental mass is strongly concentrated around the mean position, in both directions, with a larger dispersion in the streamwise direction. Experimental mass is discontinuously distributed, while the simulated results appear smoother. In this case, the mass is gradually distributed along the two directions and the maximum values are underestimated, especially for the streamwise direction.
This suggests that the model provides an acceptable simulation of the transversal diffusion, while the longitudinal one is overestimated.

Conclusions
A Eulerian model for the simulation of congested large wood transport is here proposed. It allows the one-way coupled solution of the Shallow Water Equations for the flow and the advection-diffusion equation for wood mass concentration, including proper modification to consider the floating transport. In fact, the transport velocity is expressed with an asymptotic formulation, able to mimic the inertia of large wood to the flow accelerations and decelerations, and the diffusion model comprise ad hoc diffusion coefficients.
The application of this model to a series of experiments shows a good agreement between numerical and experimental results. The model can predict the evolution of the transported wood mass and the most probable location of wood, with a percentage error of 4.5% and 16% in the streamwise and transversal directions.
Nonetheless, the diffusion appears greater than expected, with an underestimation of the maximum concentration and a consequent wider planar distribution of wood mass. This may depend on some overestimation of the diffusion coefficients or, more probably, on a numerical diffusion that remains high despite the implemented correction. An indepth analysis of the numerical diffusion tensor, comprising also different discretization methods (like finite element, e.g., [37,38]) may help assessing the actual role of the numerical diffusivity for the proposed scheme.
Overall, the model results are promising, but additional investigation is required to fix issues connected with diffusion. Additionally, since real events may include fine wood among larger wooden elements, the model behavior under such conditions should be also investigated. Finally, it would be advisable to include the effect of Large wood orientation that is a significant parameter for the application to real cases, since it affects the probability of accumulation and the consequent backwater effect.
Funding: This research received no external funding.