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

Abstract

An advection-diffusion model is proposed to simulate large wood transport during high flows. 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 flow and of the floating wood mass. The model is tested on a large series of flume experiments, involving at least 30 logs and different control parameters (flow 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 sufficient 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.

1. 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.

Focusing on wood, past and recent events demonstrated its effects both in terms of extension of flooded areas [10,11] and of potential/real damage to in-line structures, such as bridges and check dams (e.g., videos/pictures of events in Italy [12,13,14]). It is increasingly evident that wood transport should not be disregarded, especially in mountainous basins, to provide hydraulic simulations that can cope also with the residual risk connected to such phenomenon.

Recently, different models including Large Wood (LW) transport in hydraulic simulations were developed. Taking their cue from the mathematical analysis of Braudrick and Grant [15], existing models consider each log as a single element which can be entrained, transported, and deposited by the flow. These models can be defined as Eulerian-Lagrangian models, since they couple the two-dimensional solution of the Shallow Water Equations (SWE, with Eulerian methods, e.g., Finite Volumes) with a Discrete Element approach, properly Lagrangian. Wood elements are entrained and transported according to the hydrodynamic forces [16,17,18,19], computed with appropriate hydrodynamic coefficients [20,21,22], or following a kinematic approach [23]. The three-dimensional effect of the flow can also be considered, especially when focusing on the accumulation of wood at in-line structures [24].

These models were tested on different flume or field experiments [25,26,27], the most advanced ones being able to replicate the backwater effect resulting from wood clogging and the variation in flood hazard mapping [28].

Despite proving to be efficient, such models present some drawbacks. The number of logs, their initial positioning, dimensions, and shapes are required for the simulation. In addition, since several elements can be transported during a flood and each one needs to be considered in a Lagrangian framework, simulation times can grow owing to the large number of operations required to compute individual trajectories and mutual interactions. In addition, such deterministic approach is not fully representative of a stochastic phenomenon such as LW transport.

A different approach is here proposed to overcome these limits, which is particularly suitable for the congested transport condition, i.e., the transport of several LW pieces, strongly interacting and behaving like a carpet [29]. Such an approach, already drafted in [30], considers the ensemble of transported logs like a continuous quantity floating on the water surface. The total mass of the logs is thus treated as a passive and conservative substance, so that an advection-diffusion equation is numerically coupled with the SWE, allowing for the coupled solution of the water flow and the LW transport. This solution prevents the acquisition of detailed information about LW pieces that are substituted by the total mass of wood observed, or expected, during a flood.

In this contribution (i) the mathematical model outlined in [30] is discussed, and the advection-diffusion equation for wood mass transport is implemented in the two-dimensional code ORSA2D [31]; (ii) the main parameters of the advection-diffusion equation, i.e., the transport velocity and the diffusion coefficients, are calibrated to the typical situation of large wood thanks to the analysis of available experimental campaigns [30] and, finally, (iii) the performances of the proposed model are evaluated against flume experiments [30,32], providing a global measurement of how the model can replicate LW transport.

2. Materials and Methods

2.1. Mathematical and Numerical Outlines

2.1.1. 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:

$$\frac{\partial p}{\partial t}+v·\nabla p+\omega \frac{\partial p}{\partial \theta }=K_x\frac{{\partial }^{2}p}{\partial x^2}+K_y\frac{{\partial }^{2}p}{\partial y^2}+K_{\theta }\frac{{\partial }^{2}p}{\partial {\theta }^2}+S_w,$$

(1)

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:

$$p\left(\mathit{x},\theta ,t\right)=c_w\left(\mathit{x},t\right) \delta \left(\mathit{x},\theta ,t\right),$$

(2)

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:

$$\delta \frac{\partial c_w}{\partial t}+c_w\frac{\partial \delta }{\partial t}+v·\nabla (c_w\delta )+\omega \frac{\partial }{\partial \theta }\left(c_w\delta \right)=K_x\frac{{\partial }^2}{\partial x^2}\left(c_w\delta \right)+K_y\frac{{\partial }^2}{\partial y^2}\left(c_w\delta \right)+K_{\theta }\frac{{\partial }^2}{\partial {\theta }^2}\left(c_w\delta \right).$$

(3)

By integrating Equation (3) in the interval (0,2π), being ${{\displaystyle \int }}_0^{2\pi }\delta \left(\mathit{x},\theta ,t\right)d\theta =1$, due to the flow continuity equation and to the fact that c_w(x,t) and v are independent of θ, it results:

$$\frac{\partial c_w}{\partial t}+\nabla c_{w}v=K_x\frac{{\partial }^2{c}_w}{\partial x^2}+K_y\frac{{\partial }^2{c}_w}{\partial y^2},$$

(4)

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:

$$\frac{\partial \delta }{\partial t}+v·\nabla \delta +\omega \frac{\partial \delta }{\partial \theta }=K_x\frac{{\partial }^2\delta }{\partial x^2}+K_y\frac{{\partial }^2\delta }{\partial y^2}+K_{\theta }\frac{{\partial }^2\delta }{\partial {\theta }^2}.$$

(5)

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.

2.1.2. 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:

$$\frac{\partial \mathit{U}}{\partial t} + \frac{\partial \mathit{F}\left(\mathit{U}\right)}{\partial x} + \frac{\partial \mathit{G}\left(\mathit{U}\right)}{\partial y} = \mathit{S}\left(\mathit{x},\mathit{y},\mathit{U}\right),$$

(6)

being,

$$\begin{gathered} \mathit{U} ={\left(h,qx,qy,h\phi \right)}^T, \\ \mathit{F} ={\left(qx,\frac{q_x^2}{h}+\frac{g{h}^2}{2},\frac{q_x{q}_y}{h},uh\phi \right)}^T, \\ \mathit{G} ={\left(qy,\frac{q_x{q}_y}{h},\frac{q_y^2}{h}+\frac{g{h}^2}{2},vh\phi \right)}^T, \\ \mathit{S} ={\left(0,gh\left(S_{0x} - S_{fx}\right),gh\left(S_{0y} - S_{fy}\right),\overline{\nabla }\left(\mathit{K}h\overline{\nabla }\phi \right)\right)}^T, \end{gathered}$$

(7)

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₀ 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:

$$\phi =\frac{c_w{\rho }_w}{h dx dy},$$

(8)

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:

$$\begin{gathered} {\tilde{\lambda }}^1 = \tilde{\mathit{u}}·\mathit{n} + \tilde{c}, \\ {\tilde{\lambda }}^2 = \tilde{\mathit{u}}·\mathit{n}, \\ {\tilde{\lambda }}^3 = \tilde{\mathit{u}}·\mathit{n} - \tilde{\mathit{c}}, \\ {\tilde{\lambda }}^4 = \tilde{\mathit{u}}·\mathit{n}; \end{gathered}$$

(9)

$${\tilde{e}}^1 = \left(\begin{array}{c}1\\ \tilde{u}+\tilde{c}{n}_x\\ \tilde{v}+\tilde{c}{n}_y\\ \tilde{\phi }\end{array}\right), {\tilde{e}}^2 = \left(\begin{array}{c}0\\ -\tilde{c}{n}_y\\ \tilde{c}{n}_x\\ 2\tilde{\phi }\end{array}\right), {\tilde{e}}^3 = \left(\begin{array}{c}1\\ \tilde{u}-\tilde{c}{n}_x\\ \tilde{v}-\tilde{c}{n}_y\\ \tilde{\phi }\end{array}\right), {\tilde{e}}^4 = \left(\begin{array}{c}0\\ 0\\ 0\\ 1\end{array}\right).$$

(10)

In Equations (9) and (10), the eigenvalue ${\tilde{\lambda }}^4$ refers to the wood transport equation, as well as the fourth term of the eigenvectors ${\tilde{e}}^1, {\tilde{e}}^2, {\tilde{e}}^3$ and the eigenvector ${\tilde{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 φ:

$$\begin{gathered} \tilde{\mathit{u}} = \frac{u_R\sqrt{h_R}+u_L\sqrt{h_L}}{\sqrt{h_R}+\sqrt{h_L}} , \tilde{v} = \frac{v_R\sqrt{h_R}+v_L\sqrt{h_L}}{\sqrt{h_R}+\sqrt{h_L}}, \\ \tilde{\mathit{c}} = \sqrt{g\frac{h_R+h_L}{2}} , \tilde{\phi } = \frac{{\phi }_R\sqrt{h_R}+{\phi }_L\sqrt{h_L}}{\sqrt{h_R}+\sqrt{h_L}}, \end{gathered}$$

(11)

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).

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:

$${{\displaystyle \iint }}_{{\Omega }_i}^{ }\overline{\nabla }\left(\mathit{K}h\overline{\nabla }\phi \right) d\Omega = {{\displaystyle \oint }}_{\partial {\Omega }_i}^{ }\left(\mathit{K}h\overline{\nabla }\phi \right) \mathit{n} ds = {\displaystyle \sum }_{k=1}^{NE}{\left(\frac{{\mathit{K}}_R+{\mathit{K}}_L}{2}\mathit{n} min\left(h_R^{\star },h_L^{\star }\right)\frac{\delta \phi }{d_{RL}}\right)}_k{\mathit{n}}_k{s}_k,$$

(12)

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:

$$K_{corr}=\frac{1}{2} \left|\mathit{u}·\mathit{n}\right| \left(d_{RL}-\left|\mathit{u}·\mathit{n}\right|\Delta t\right),$$

(13)

and it is implemented as follows:

$${\displaystyle \sum }_{k=1}^{NE}{\left(K min\left(h_R^{\star },h_L^{\star }\right)\frac{\delta \mathit{\varphi }}{d_{RL}}\right)}_k{s}_k,$$

(14)

with

$$\{\begin{array}{cc}K=0& if K_{corr}>{\left(\frac{{\mathit{K}}_R+{\mathit{K}}_L}{2}\right)}_k{\mathit{n}}_k{\mathit{n}}_k\\ K=max\left({\left(\frac{{\mathit{K}}_R+{\mathit{K}}_L}{2}\right)}_k{\mathit{n}}_k{\mathit{n}}_k-K_{corr},0\right)& if K_{corr}\le {\left(\frac{{\mathit{K}}_R+{\mathit{K}}_L}{2}\right)}_k{\mathit{n}}_k{\mathit{n}}_k\end{array}$$

(15)

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)).

$$K_s = 0.005 {\left(F{r}_p\right)}^{-0.740}{\left(t_{Rr}\right)}^{-0.300}, K_t= 0.011 {\left(F{r}_p\right)}^{0.188}{\left(t_{Rr}\right)}^{0.271};$$

(16)

$$K_{xx}=\frac{K_{s}u+K_{t}v}{\sqrt{u^2+v^2}},K_{yy}=\frac{K_{s}v+K_{t}u}{\sqrt{u^2+v^2}},K_{xy}=\frac{1}{2}\left(K_s-K_t\right)\frac{2uv}{u^2+v^2};$$

(17)

where $F{r}_p = \frac{v_{MAX}}{\sqrt{g{L}_w}}$, $t_{Rr} = \frac{t_R}{L_W}$, 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.

2.2. 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₅₀ = 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.

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⁻³, such as that of green wood broadleaved trees [36].

Table 1. Summary of the main parameters of the experiments.

Flume	Slope i [−]	Discharge Q [L s−1]	Froude Number Fr [−]	Release Distance tR [m]	Log Length LW [m]	Log Diameter DW [m]	Number of Repetitions NL
S-shaped	0.0004	20	0.1	0.8	0.3	0.02	41
					0.4		30
		50	0.2	0.4	0.3		35
					0.4		46
				0.8	0.3		29
					0.4		41
	0.0016	81	0.4	0.8	0.2	0.014	40
						0.02	46
					0.3	0.02	42
					0.4	0.02	40
						0.03	39
S-shaped with narrowing	0.0016	66	- 1	0.4	0.2	0.014	46
						0.02	42
					0.3	0.02	48
					0.4	0.02	47
						0.03	49
				0.8	0.2	0.014	44
						0.02	42
					0.3	0.02	48
					0.4	0.02	48
						0.03	45
				1	0.2	0.014	48
						0.02	48
					0.3	0.02	54
					0.4	0.02	48
						0.03	49

¹ Variable along the flume.

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=\frac{v_{mean}}{\sqrt{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.

2.3. 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:

$$v_{MAX} = \frac{1}{0.85} v_{mean}.$$

(18)

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⁻¹ 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:

$$c_{vel} = \left(0.3+0.19\left(1-e^{-\frac{t}{c}}\right)\right)\frac{1.18}{0.49},$$

(19)

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:

$$c=4.39F{r}_R.$$

(20)

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:

$$k_d=1+0.4t_R.$$

(21)

Overall, the transport velocity computed at each time step, for each cell i containing a wood concentration, is:

$$\{\begin{array}{cc}v\left(i,t\right)=\left(0.3+0.19\left(1-e^{-\frac{t}{c}}\right)\right)\frac{1.18}{0.49} & if v_f\left(i_{c max},t\right)-v_f\left(i_{c max},t-1\right)>\Delta v\\ v\left(i,t\right)=\left(0.3+0.19\left(1-e^{-\frac{t}{c}}\right)\right)\frac{1.18}{0.49}{k}_d& if v_f\left(i_{c max},t\right)-v_f\left(i_{c max},t-1\right)\le \Delta v,\end{array}$$

(22)

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⁻¹ 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).

2.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:

$$\frac{\partial \phi \left(\mathit{x},t\right)}{\partial t}=\frac{{\rho }_w{{\displaystyle \sum }}_{i=1}^{N_L}{L}_w{}_i\frac{D_w{}_i^2}{4}\pi }{A_{cell}{h}_{cell}{T}_{ref}}$$

(23)

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].

3. Results

3.1. 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² = 0.71) than for u component (Figure 5a, average correlation coefficient R² = 0.20), which is less variable in the simulation than in the experiments (especially for the lower Froude numbers).

Figure 6a shows the computed and measured Froude number along the S-shaped flume with narrowing, while Figure 6b,c show the comparison of the flow velocity. The correlation coefficients are 0.79 for the Froude number, 0.77 for the u component and 0.69 for the v component, confirming the good representation of the flow in this flume. Most likely, the fact that the downstream boundary condition was based on experimental measures affected the higher quality of the numerical simulation, compared with the test in the S-shaped flume.

3.2. 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.

3.2.1. S-Shaped Flume

Graphs in Figure 7 show the comparison of the experimental log position and the computed wood concentration for three tests in the S-shaped flume, with different flow conditions. In each sub-figure, a time sequence is shown, with a time interval of about 20 s, 10 s and 5 s, shorter for the test with higher Froude number. The minimum value of the volume averaged concentration that corresponds to the largest contour line, is 5 kg m⁻³ and the last time instant varies depending on the considered test (140 s, 91 s and 75 s for the shown graphs).

3.2.2. 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⁻³. 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).

4. Discussion

Figure 7 and Figure 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.

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 Figure 7 and Figure 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⁻³) 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 Figure 7 and Figure 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. Figure 9a and Figure 10a show that the mean planar position for each test is, on average, in the range 14.5–16.5 m in the streamwise direction for the S-shaped flume and 15.5–16.5 m for the flume with narrowing (with one exception for both configurations). Average transversal positions are grouped around the flume axis (+0.2 m, −0.6 m). The numerical maximum, observed at the same time selected for the experiments, is scattered longitudinally for the S-shaped flume (13–17.5 m), while it is well grouped for the flume with narrowing, although in a lower range of streamwise coordinates (14.75–16 m). In both cases, the transversal ranges are like the experimental ones.

The differences found for the two flume configurations are confirmed by Figure 9b and Figure 10b, which show the distance vector between the numerical and experimental average positions, giving a visual representation of the error in estimating the most probable position of the transported logs. For the S-shaped flume (Figure 9b), in most cases these vectors point in the downstream direction, showing that the numerical concentration tends to flow faster than experimental logs. On the contrary, the numerical concentration is slower in the flume with narrowing, with vectors pointing upstream (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\raisebox{1ex}{$D_W{}^2$}\!\left/ \!\raisebox{-1ex}{$4$}\right.{\rho }_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 Figure 7 and Figure 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.

Figure 12 and Figure 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 Figure 12 and Figure 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⁻².

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.

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.

5. 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 in-depth 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.