Two‐Dimensional Numerical Modeling of Large Wood Transport in Bended Channels Considering Secondary Current Effects

The modeling of large wood (LW) transport in rivers has received increasing interest from researchers in the last decade due to the widely recognized role of LW concerning flooding risk. For this purpose, few 2D depth‐averaged hydraulic models have been coupled with LW transport models. However, such models usually neglect the effects of secondary currents on LW trajectories in river bends. In this work, the model Iber‐Wood was enhanced to simulate the effects of secondary currents in river bends on LW trajectories. The proposed methodology presents a new formulation for considering secondary current effects on the flow field derived from the Manning formula and considers a new approach for reproducing the surface flow field that develops at channel bends. The enhanced model was tested to reproduce a series of laboratory experiments on wood transport in a sharp channel bend. The methodology introduces two new parameters in the model related to the secondary current effects, that is, the secondary current intensity and the adaptation length. These parameters were calibrated using available data from laboratory experiments. The good agreement between observed and simulated dowel trajectories in a sharp channel bend validated the proposed approach to simulate LW transport in the case of secondary currents.

Recently, numerical modeling tools have been developed to explicitly simulate wood transport at the river reach-scale (Bertoldi & Ruiz-Villanueva, 2017;Kang & Kimura, 2018;Mazzorana et al., 2011Mazzorana et al., , 2018;;Persi et al., 2018;Ruiz-Villanueva, Bladé, et al., 2014;Ruiz-Villanueva et al., 2020;Zischg et al., 2018).These models generally solve the two-dimensional Shallow Water Equations (2D-SWE) for modeling hydrodynamics and use different approaches for sediment and wood transport.However, until now, only two existing models in the literature fully couple wood transport to hydrodynamics (Kang & Kimura, 2018;Ruiz-Villanueva, Bladé, et al., 2014).One of the first attempts to simulate wood transport in rivers was the model presented by Mazzorana et al. (2011) which considered the whole process chain of LW recruitment, transport, and deposition.The model enabled the calculation of LW pathways for a given wood volume and computed the transport conditions using results from a 2D hydrodynamic model.In this case, the one-way coupling model considered the flow forces acting on wood elements to calculate LW movement, but the influence of wood on the hydrodynamics was neglected.Similarly, Zischg et al. (2018) presented a modified version of Mazzorana's one-way coupling model for simulating the recruitment, transport, and deposition of LW at a specific point in a river system.In this approach, the influence of wood on the hydrodynamics remained neglected.
The interaction between water flow, wood, and morphodynamic processes was implemented by Ruiz- Villanueva, Bladé, et al. (2014) and Kang and Kimura (2018).Both 2D models proposed a two-way coupling approach between a 2D-SWE Eulerian model, used for hydrodynamics, and a 2D Lagrangian or discrete element model for simulating the motion of individual pieces of wood.Kang and Kimura (2018) presented a two-way coupling model in which the mutual interaction between LW and the flow was modeled by calculating the drag force exerted by the wood on the flow.This methodology was tested by simulating flume experiment observations on the transport and deposition of LW pieces with and without roots.The two-way coupling model Iber-Wood presented by Ruiz- Villanueva, Bladé, et al. (2014) was also based on an additional drag force.Iber-Wood has been extensively tested, validated, and applied to study various aspects related to wood dynamics in rivers with different morphologies (e.g., Ruiz Villanueva, Bladé Castellet, et al., 2014;Ruiz-Villanueva et al., 2016, 2017).However, the model has been less often applied to study LW transport in meandering or bended channels.In such morphologies, wood transport might be significantly affected by the secondary currents that appear in bends, which may influence the wood trajectories, especially when floating.This important aspect has been largely overlooked until today.This work aims to fill this gap by adapting the numerical model Iber-Wood to simulate secondary currents and, thus, to consider the effect on the LW transport.For the validation of the enhanced model, we used available data from flume experiments.

Large Wood Transport in Channels With Bends: Theoretical Background
The curvature of river bends induces an imbalance of the local pressure gradient and the centrifugal force followed by the establishment of a secondary current, also called spiral flow or helical flow (Rozovskii, 1957).This secondary flow is associated with the streamline curvature and is typically directed toward the outer bank at the free surface.In meandering channels, a further mechanism needs to be considered due to the topographically driven secondary flow (Dietrich & Smith, 1983).The topographic steering of the flow arises in case of longitudinal variations of bed elevation and longitudinal velocity which, by simple continuity requirements, produce an additional lateral velocity.
The streamline curvature results in a three-dimensional flow that induces a redistribution of velocities, a deviation of boundary shear stress directions, and causes additional friction losses (Blanckaert & De Vriend, 2003).The topographic steering, on the other hand, results in a mass redistribution that determines the flow of mass movement over the deepest zones of the bend (Blanckaert, 2010).Globally, at the surface, where LW is commonly transported as floating objects (Diehl, 1997;Lagasse et al., 2010), the boundary shear stress is directed toward the outer bank of the channel and contributes, together with the effect of the wood piece inertia (Persi et al., 2018), in generating the driftwood (Diehl, 1997;Schmocker & Weitbrecht, 2013).This mechanism, depicted in Figure 1, drives the accumulation of LW typically at the outer, downstream bank of channel bends or meanders (Abbe & Montgomery, 2003;Gurnell et al., 2002;Le Lay et al., 2013).

Numerical Modeling of Secondary Currents
Due to the three dimensionality of secondary current flows in natural rivers, 3D models are required to properly simulate secondary currents.However, applying 3D models to natural river domains is still a very time-and computational-consuming procedure (Deltares, 2019;Lane et al., 1999;Nabi et al., 2016), and only a few attempts have been made to apply 3D models to simulate LW transport (e.g., Kimura et al., 2021).2D-SWE-based models can be adapted to simulate the effects of the secondary current on the mean flow, that is, to consider the additional stresses due to helical flows on the depth-averaged velocity field, or to estimate the surface or near-bed velocities, with a good compromise between accuracy and feasibility (Baghlani, 2012;Ghamry & Steffler, 2005;Nabi et al., 2016;Rinaldi et al., 2008;Vasquez et al., 2005).
The integration of the three-dimensional continuity equation and momentum equations (Reynolds Averaged Navier-Stokes [RANS] equations) from the bed to the water surface provides the 2D-SWE.The resulting equations can be expressed, considering a Cartesian frame of reference, as (Blanckaert & De Vriend, 2003) where t is time (s), h is the water depth (m), U x and U y are the depth-averaged velocities in x-and y-directions (m s −1 ), respectively, g is the gravitational acceleration (m s −2 ), z b is the bed elevation (m), τ bx and τ by are the shear stress components at the bed in x-and y-directions (kg s −2 m −1 ), respectively, ρ is the water density (kg m −3 ), and F sx and F sy are dispersion terms (m s −2 ).
The dispersion terms result from the vertical nonuniformity of velocity distributions.They are usually negligible in nearly straight reaches (F sx ≈ F sy ≈ 0) but are of importance in bends due to the formation of secondary currents (Baghlani, 2012;Deltares, 2019;Nabi et al., 2016).The evaluation of these dispersion terms, therefore, represents the major challenge for considering the effects of the secondary current on the mean flow.In the past, few researchers provided analytical approaches to nearly horizontal flow, such as Kalkwijk and Booij (1986) and Odgaard (1986).The approach presented by Kalkwijk and Booij (1986) is based on the analysis of the convection of momentum of secondary currents in streamwise direction.It considers both the effects of Coriolis acceleration and channel curvature.This approach was applied in recent works (Baghlani, 2012;Nabi et al., 2016;Rinaldi et al., 2008).Kalkwijk and Booij (1986) considered the streamwise lag-effect on the secondary current by introducing the concept of "adaptation length" or "relaxation length," that is, the length after which the flow adapts to the change in curvature of the channel, to represent the grow/decay factor of the secondary current.Figure 1 shows the adaptation length (λ [m]) for the growth of a fully developed helical flow indicated as "generation length."In contrast, the "relaxation length" indicates the distance after which the helical flow vanishes downstream of the bend.Jagers (2003) proposed a depth-averaged advection equation for simulating the streamwise lag-effect in 2D-SWE models based on the Kalkwijk and Booij's approach.The equation proposed by Jagers (2003) allowed the author to calculate the variation in space and time of the secondary current intensity, that is, the magnitude of the secondary current component normal to the depth-averaged flow, as where I is the secondary current intensity (m s −1 ), and S I is the source term of the secondary current intensity (m s −2 ); further details regarding this equation are reported in Section 2. When the secondary current intensity is known, it can be used for considering the effects on the momentum equations and for computing bed load transport.
All the above-mentioned works aimed at simulating the effects of secondary currents on the main flow in terms of water depths and velocities (Baghlani, 2012;Deltares, 2019;Jagers, 2003;Kalkwijk & Booij, 1986;Nabi et al., 2016) and on bedload transport (Rinaldi et al., 2008;Vasquez et al., 2005).However, according to the author's knowledge, to date, there are no studies on the use of 2D-SWE models in simulating the effects of secondary currents on the transport of LW.

Iber-Wood Model
Iber-Wood is a two-way model that couples the hydrodynamic Eulerian model Iber (Bladé et al., 2014) to a LW Lagrangian or discrete element model (Ruiz- Villanueva, Bladé, et al., 2014).Iber computes the hydrodynamics by solving the 2D-SWE with or without considering turbulence and sediment transport (see Cea at el., 2007, for more details).It uses the finite volume method with a time-explicit second-order and nonoscillatory extension of Roe's upwind scheme on nonstructured meshes (see Bladé et al., 2014 for details).Iber-Wood simulates the transport or advection of individual LW elements by fully coupling the LW and hydrodynamic models: the hydrodynamic solution is used by the Lagrangian LW model to simulate the movement of the wood element that introduces additional friction and associated energy dissipation into the 2D-SWE.The LW incipient motion is based on a balance of forces acting on the LW center of mass.The LW transport simulation considers two possible mechanisms based on wood density and flow conditions: (a) floating at the water surface or (b) sliding or dragging on the riverbed.In the case of floating conditions, the transport simulation is carried out by performing a kinematic analysis, that is, the LW velocity is assumed to be the same as the flow velocity unless turbulence is computed.On the contrary, in the case of sliding conditions, the friction force is the main parameter driving the movement of the element, and the LW velocity is different from the flow velocity.In addition to translation, the model considers LW rotations and LW interactions with other wood elements or river boundaries (i.e., riverbed morphologies and in-channel structures).Further details are provided by Ruiz- Villanueva et al. (2020).
Iber-Wood allows to define initial and boundary conditions for LW.The use of initial conditions permits the user to define the entry time step for any single wood element, its position and orientation, its dimensions, and its density.On the other hand, boundary conditions can be applied only to the computational domain boundaries as a LW rate per minute.In the latter case, the wood element characteristics (dimensions and density) and its initial orientation result from a stochastic choice in the range of values defined by the user.

Implementation of River Bend Secondary Current Effects in Iber-Wood
The adaptation of Iber-Wood for simulating the effects of secondary currents due to the curvature of river bends on the LW transport allows the model to consider (a) the effects on the mean flow, (b) the effects of the classical redistribution of velocities at the water surface, and (c) the streamwise lag-effect on the secondary current.The implemented methodology was developed for a computational domain discretized with a mesh of triangles, quadrilaterals, or its combination, as Iber-Wood is based on the finite volume method.
A total of three methods were implemented in the model for this study: the new proposed methodology in the following M1; two methods selected from the literature (Kalkwijk & Booij, 1986;Odgaard, 1986) in the following M2 and M3, respectively.M2 and M3 were considered in the present study as reference methods and used as a term of comparison for measuring the performances of the proposed methodology (M1).
The main characteristics of the methodology used for adapting Iber-Wood to simulate the effects of secondary currents on the flow and, therefore, on the LW transport are summarized in Figure 2.

Proposed Methodology
In the new proposed methodology (M1), the horizontal momentum equations were modified by introducing the dispersion terms F sx and F sy as shown in Equations 2 and 3.These terms were derived by the Manning equation (Manning, 1891), making some assumptions: • the effects of the secondary current near the riverbed are comparable to those at the free water surface, so the velocity field at these two boundaries can be considered equal in terms of absolute velocities (i.e.,  ⃖⃖⃖ ⃗ b = ⃖⃖⃖ ⃗ s , where b and s subscripts mean near-bed and surface, respectively), while the direction of velocity vectors is the opposite; • at the water surface, the secondary current deviates the velocity vectors by an angle δ with respect to the mean flow velocity field, as • the surface velocity field differs from the mean flow velocity field due to the presence of the secondary current, thus, the velocity vector can be expressed as where  ⃖⃖⃖ ⃗  ′ is the velocity vector component due to the secondary current (m s −1 ); • the hydraulic radius of the cells of a 2D mesh can be approximated by the value of the water depth in the cell.
Considering these assumptions, the Manning equation can be written for the x-and y-directions as 10.1029/2022WR034363 6 of 16 where S x and S y are the slopes of the hydraulic grade in x-and y-direction (-), respectively, and n is the roughness Manning coefficient (s m −1/3 ).The slopes of the hydraulic grade are defined by two factors: one is related to the secondary current (i.e.,   ′  and   ′  in Equations 8 and 9, respectively).The other refers to the mean flow velocity.The implemented method for calculating F sx and F sy considers only the term related to the secondary current, resulting in in which the velocity vector components due to the secondary current (   ′  and   ′  ) are calculated from Equation 7 according to The streamwise lag-effect on the secondary current is modeled by using Equation 4. This equation requires the definition of the secondary current intensity I and the source term S I .The intensity I changes across the channel bend depending on the local radius of the curvature of the streamline, water depth, and flow velocity, and its equilibrium value can be expressed as (Deltares, 2019;Jagers, 2003;Nabi et al., 2016) where  ⃖⃖ ⃗  is the depth-averaged velocity vector (m s −1 ) and R s is the radius of curvature of the streamline (m) that is considered as the curvature of flow streamlines by means of the velocity field.It is defined as The source term of the secondary current intensity is defined according to Jagers (2003) as where I FD is the spiral flow intensity for a fully developed constant radius bend (m s −1 ) and λ is the adaptation length (m) that can be computed using the formulation proposed by Kalkwijk and Booij (1986).
where α is a nondimensional parameter that depends on the von Karman constant κ (-) and the Manning coefficient.Alternatively, it can be set by the user.In this case, the user can specify the generation length (λ gen ) and the relaxation length (λ rel ) that control respectively the generation and the decay of spiral flow according to The solution of Equation 4 provides the variation of the secondary current intensity in space and time.Next, the value of I is used to determine the velocity components at the water surface from Equations 5 and 6.The angle δ is defined as where c is a nondimensional term (-), A δ is a nondimensional parameter that depends on the von Karman constant κ (-) and the Chezy coefficient 2 ), and   * s is the effective radius of curvature (m) calculated from Equation 14 by knowing the secondary current intensity (   * s = ℎ| ⃖⃖ ⃗  |∕ ).The correction term c is inserted in the definition of δ to account for the spiral flow characteristics.The shear stress deviation is not constant along the cross section since it is practically zero close to riverbanks (where the flow is mainly vertical) and increases moving from the riverbanks to the main active channel (Blanckaert & De Vriend, 2003).To reproduce this effect, the correction term is set equal to zero next to the banks or channel walls and calculated elsewhere as where c m is defined for each mesh element at each time step as the minimum value of the correction term c of border mesh elements (-), d is the distance to the boundary or lateral channel wall (m), and ∆c is the gradient of the correction term (-).The gradient ∆c expresses the variation of c along the channel cross section.It is set by the user and maintained fixed during the simulation.
This methodology introduces new terms in the model related to the simulation of the secondary current effects, among which the most important are λ and ∆c.The adaptation length ∆c can be calculated with Equation 17or defined by the user (Equation 18), while the gradient of the correction term ∆c always needs to be defined by the user.Both terms, or only ∆c, can be used by the user as calibration parameters for the simulation of secondary current effects.

Reference Methodologies
Methods M2 and M3 are briefly introduced here, and further details are given in Supporting Information S1.The method M2 follows the adaptation of secondary currents in nearly horizontal flow proposed by Kalkwijk and Booij (1986).The main flow is assumed to have a logarithmic velocity profile, and the secondary current originates from multiplication of a universal function with the spiral motion intensity.The method M3 is based on the analytical approach presented by Odgaard (1986), who presented a model for steady, subcritical, turbulent flow in alluvial channel bends with uniform bed sediment.A definition of the additional dispersion terms F sx and F sy (see Equations 2 and 3) can be derived from the Odgaard model.For both M2 and M3, the streamwise lag-effect on the secondary current was modeled by using Equation 4.

Experimental and Numerical Setup
To study the trajectories of single wooden dowels in sharp bends, flume experiments were carried out in a double-bended channel in the laboratory of the Leichtweiß-Institute for Hydraulic Engineering and Water Resources of the Technische Universität Braunschweig, Germany.These experiments, detailed in a previous publication by Innocenti et al. (2022), involved varying wood dimensions (i.e., diameter and length), wood inlet points, and initial wood orientations.Here, we provide a brief summary of key aspects, with further details available in the referenced publication.
A channel 26 m long, 2.40 m wide, and 0.40 m deep had a rectangular cross-sectional channel, and two bends were located 5.46 and 13.20 m downstream from the inlet.Laboratory observations were limited to the first bend, a bend to the left characterized by a curvature ratio of 1.5.Experiments were performed by reproducing the flow conditions previously studied by Zaid (2018) who, using the same channel, characterized the hydraulics at the first bend of the channel by measuring superelevations and the flow field velocity.These experiments were carried out with a constant discharge of 0.13 m 3 s −1 and a water depth of 0.10 m along the channel.The resulting Froude number was 0.55.Velocity measurements were performed at nine cross sections along the bend (see the top view in Figure 3).For each cross section, seven vertical profiles were measured consisting of nine measurement points (from 0.01 to 0.09 m above the bottom with a spacing of 0.01 m, see the cross section view in Figure 3).The velocity measurements were also used to calculate the flow streamlines.Water depths were measured using point gauges located at cross sections 1, 5, and 9 and in the straight reaches (see the top view in Figure 3).
Cylindrical wooden dowels (i.e., branches and roots were neglected) with four different sizes were used in the experiments to simulate LW and explore the effect of the size on the simulation of LW trajectory.The present study specifically focuses on the observation of the 0.40 m long dowel with a diameter of 0.06 m, dowel type IV in Innocenti et al. (2022).The dowels were inserted individually by hand upstream to the bend considering three inlet points (P1, P2, and P3, respectively; see the top view in Figure 3).Additionally, two initial dowel orientations were considered: parallel and perpendicular to the channel axis.In summary, data of six different experiments were available in which the dowels were tracked to obtain their trajectories.For further details on dowel-experiments, refer to Innocenti et al. (2022).
For the numerical simulations, the channel geometry was discretized using a structured mesh with a mean cell size of 0.1 × 0.1 m based on the channel and dowel dimensions.The hydraulic model boundary conditions were set according to laboratory experiments: a constant inlet discharge of 0.13 m 3 s −1 was defined at the upstream boundary, while a fixed water depth of 0.10 m was set at the downstream boundary.LW elements were entered into the model by defining a wood initial condition reproducing the laboratory experimental setup.The turbulence was modeled using the k-epsilon model (Rastogi & Rodi, 1978) already implemented in Iber.Calibration of hydraulics, including the Manning coefficient and the new secondary currents terms ∆c and λ, was performed by comparing numerical results with all available observations of water depth and surface velocity before performing dowel simulations.The final values of the calibrated parameters were n = 0.01 s m −1/3 , ∆c = 0.4, and λ = 3 m for both generation and relaxation lengths.

Model Performance
The experiments were simulated with all three methods implemented and using Iber-Wood without any correction for secondary effects of the current.The model performance was evaluated by comparing the numerical results with laboratory observations (i.e., water depth, flow velocity, and dowel trajectories).The ability of the methods to reproduce the observed hydraulic characteristics was tested by calculating the root mean square error (RMSE), the relative standard error (RSE), and the Nash-Sutcliffe efficiency (NSE) coefficient (Nash & Sutcliffe, 1970) and by measuring the similarity of the surface flow lines by defining a coefficient of similarity (C s ).The latter was also used to assess the similarity between observed and simulated wood trajectories.The 10.1029/2022WR034363 9 of 16 similarity coefficient (C s ) is a dimensionless coefficient that expresses the similarity of the observed and simulated flow lines or LW trajectories in a predefined domain (i.e., in a part of the channel bounded by reference cross sections).C s was calculated using the ratio between the absolute value of the horizontal area between the observed and simulated streamlines or trajectories and the total area of the predefined domain (A tot [m 2 ]) according to the following formula: where A i,obs and A i,sim denote the observed and simulated horizontal areas of the pin lines or trajectories (m 2 ).A i,obs and A i,sim were calculated as the areas under the curves (streamlines or trajectories) in a horizontal plane.that is, in the local coordinate system.The value of A tot was set to be equal to the area of the channel between cross sections 1 and 9 when calculating the similarity of the streamlines and to be equal to the channel area between cross sections CSS and CSE when calculating the similarity of the LW trajectories (see Figure 3 for the location of reference cross sections).
In general, the coefficient C s ranges between 0 (i.e., maximum similarity, the two lines overlap) and 1 (i.e., worst fit).
In addition, the ability of the model to properly reproduce the observed dowel trajectories was evaluated by considering the travel time.For each dowel, the travel time was calculated as the required time to float from cross section CSS to CSE.

Modeling the Effects of the Helical Flow
Figure 4 reports the results of the simulation in terms of the secondary current intensity (I) and the deviation angle (δ) values using method M1.The intensity (Figure 4a) resulted in a central region in which the effects of secondary current were higher and varied along both streamwise and transverse directions.In the streamwise direction, the change in curvature of the channel induced the generation of the secondary current that is reflected by the increased intensity.At the end of the bend, I progressively decreased once the curvature changed again.Thus, the presence of the secondary current (represented through its intensity) was not only limited to the bend but also affected areas closer upstream and downstream of the bend.In the transverse direction, I progressively decreased moving from the inner to the outer parts of the channel.Consequently, the deviation angle at the surface (Figure 5b) changed accordingly to the variation of I.The angle δ was maximum at the central region of the channel bend that extended approx imately from cross section 2 to cross section 8.In cross-sectional direction, the mesh cells close to the channel walls were characterized by an δ equal to zero that prevented vectors from being directed outward from the channel.

Effects on the Hydrodynamics
The comparison between observed and simulated values of water depth and water surface velocity magnitude is shown in Figure 5 and the performance of the considered methods is reported in Table 1.Values for the individual cross sections are reported in Supporting Information S1.
Overall, comparable performances were observed for the model without any correction and methods M1 and M3, for which the NSE ranged between 0.60 and 0.65, and between 0.26 and 0.33, in reproducing water depth and surface velocity observations, respectively.Method M2 produced slightly different results, particularly in simulating the water depth (NSE = 0.1), while it provided the best agreement with surface velocity observations (NSE = 0.41).The difference between M2 and the other methods could be noted also from the RMSE and the RSE which were observed to be about 75% and 15% greater for M2 in reproducing the water depth and the surface velocity, respectively.However, RSE values were in the range of 5%.
Figure 6 shows the comparison of the surface flow streamlines and can be used to investigate the effects of the implemented methodologies on the surface velocity flow field.The streamlines s1-s3 were derived considering three starting points at cross section 1.The streamline starting points at cross section 1 were selected to be representative of the lateral and central areas of the channel.Qualitatively, the differences between experimental observations and models increased in the second part of the bend (i.e., between cross sections 5 and 9).These differences were quantified by calculating C s by using the observed streamlines as the target.Figure 6b shows the measured C s values for the entire length of the bend (i.e., from cross sections 1 to 9).The proposed method M1, for which the lowest C s values were calculated, resulted in the best match.On the other hand, the model without any corrections for secondary current effects was characterized by the highest C s values.
The similarities for M2 and M3 were in between the other two methods (the results suggest that M2 was closer to the observation than M3).Moreover, the range of variation in C s between the four methods was the lowest for the streamline s3 and the highest for the streamline s2.
Figures 6c and 6d report the measured similarity for the first (i.e., from cross sections 1 to 5) and second (i.e., from cross sections 5 to 9) parts of the bend, respectively.Considering the two parts of the bend separately, the results confirmed what was observed qualitatively for the entire bend.In the first part (Figure 6c), the range of variation in C s between the four methods was lower than in the second part (Figure 6d).In addition, in the first part, the performances for methods M1 and M2 were comparable, with M2 showing the best similarity for the streamline at the center of the channel (s2).In the second part of the bend, the similarities decreased for the four methods, resulting in the highest C s values.No differences were observed in this second part from what was observed for the entire bend.

Large Wood Trajectories With Secondary Currents
Figure 7 shows the LW trajectory results of the dowel-experiments simulated with Iber-Wood.Figures 7a and 7b report the results for the experimental configuration in which the dowel was inserted from P 2 with the dowel initial orientation parallel to the channel axis.The similarity C s for the six considered experiments and the four models is reported in Figure 7c.The observed trajectory (see Figure 7a) drifted toward the out wall of the channel once the dowel entered the bend.This mechanism, that is, the dowel drift, was not reproduced by the model without any corrections for which the similarity coefficient (C s ) resulted in the highest for the four considered methods.On the contrary, the three implemented methods allowed the reproduction of the dowel drift, bringing the simulated trajectories closer to the observed one.As for the streamline comparison (see previous section), the methods M1 and M2 performed similarly, although M1 performed slightly better.The performance of M3 was in between M1/M2 and the model without corrections.Method M1 performed best for all the considered experimental configurations from inlets P 1 and P 2 , while this was not the case for inlet P 3 (see Figure 7c).In general, for the configurations where dowels were released from the inlet P 1 , the similarity was less than for the other inlet points and initial dowel orientations considering the four used methods.Contrarily, the similarity of all methods became comparable for inlet point P 3 .In fact, considering P 3 the range of variation in C s between the best and the worst method was reduced by about 6 times with respect to what was measured for P 1 and P 2 .In addition, the results of the model without any correction were comparable to the ones obtained with implemented methodologies only for P 3 experiments, while they resulted the worst for all the other experiments.
Figure 8 shows the comparison between observed and simulated dowel travel times from the cross section CSS to cross section CSE (see Figure 3) showing that the numerical simulations underestimated the dowel travel times.For all configurations, the inlet point strongly determines the travel time of dowels both in experiments and simulations, resulting in the longest travel time for P 3 and the lowest for P 1 .The best correlation was obtained using the method M1 and the worst using the model without corrections.

Sensitivity Analysis
The sensitivity analysis of the proposed methodology (M1) to the variables related to the secondary current effects (i.e., ∆c and λ) is reported in Figure 9.The method appeared more sensitive to the gradient of the correction term   ∆c (Figures 9a and 9b) than to the adaptation length λ (Figures 9c and 9d).The simulated trajectories obtained varying λ did not change much (i.e., C s ranged between 0.011 and 0.013), either by considering Equation 17to calculate the variable or by imposing a fixed value for λ gen and λ rel .Moreover, not even an increase of an order of magnitude induced significant differences, resulting in a C s variation of 0.002 while using λ gen and λ rel equal to 1 or 10.Contrarily, the change of ∆c produced significant effects on the dowel trajectory, that is, increasing the dowel drift effect moving downstream along the bend.In this case, increasing ∆c from 0.1 to 0.8 produced an increment of the similarity of about 77%.

Methodology Implications and Limitations
The implemented methodologies for the simulation of secondary currents at the water surface were found to be crucial for accurately predicting the LW trajectories along river bends.This outcome is based on the comparison of observed and simulated trajectories, as well as the analysis of dowel travel times (see Figures 7 and 8).However, some limitations of the proposed methodology deserve attention.
A difference was observed considering experiments from the inlet point P 3 .As highlighted by Innocenti et al. (2022), the dowels inserted from the inlet point P 3 frequently interacted with the right wall of the channel while moving downstream along the bend.These interactions, which are important in determining the dowel trajectory in the laboratory, were not reproduced by any of the considered methods.The discrepancy arises from the definition of the correction term c that appears in calculating the deviation angle at the surface δ (see Equation 19).The term c was set equal to zero next to the channel walls, where the helical flow is predominantly vertical.While this assumption prevents outward-directed velocity vectors, it also results in wood trajectories parallel to the channel wall, thus not allowing interactions.
The gradient ∆c was calibrated in this study using available hydraulic data of flow depth and velocity.This approach is convenient for cases where flow measurements are available, otherwise, the definition of this parameter must be made based on literature data and, in any case, will depend on the user's experience.For the specific case of a sharp river bend (curvature ratio less than 2), as defined by Blanckaert et al. (2013), ∆c = 0.4 is recommended.
In addition, the study reports tests conducted considering a channel with a rectangular cross-sectional geometry and fixed bed conditions.For these reasons, the role of the topographic steering on the LW transport prediction was completely neglected.It is possible to assume that in the case of natural bends or meanders, the cross-sectional shape of the river induces a change in the LW transport dynamics due to the different velocity redistribution.Accordingly, modifications to the presented methodology, particularly regarding the definition of the deviation angle δ, will be necessary to adequately reproduce the surface flow field while considering natural bends or meanders (Blanckaert, 2010;Diehl, 1997;Innocenti et al., 2022;Panici, 2021).

Conclusions
The study aimed to numerically reproduce the effects of secondary currents on wood transport in the existing 2D model Iber-Wood.The model underwent enhancement through the implementation of a methodology designed to simulate the effects of the helical flow on the surface flow field.The effects on the flow field were addressed through the definition of a methodology derived from the Manning formula representing a novel approach to replicating the surface flow field that develops at channel bends.The performance of the proposed method was tested by reproducing laboratory experiments.The resulting agreement between observed and simulated dowel trajectories in a sharp channel bend confirmed the reliability of the proposed approach.
At present, this methodology stands as the first attempt to incorporate secondary current effects in a 2D depth-averaged model for simulating the transport of LW in bended channels.Further application of the method, including exploration of different channel geometries or its extension to river-scale applications, will facilitate refinement and enhancement of the methodology.This promises to provide a robust tool for simulating LW transport in rivers for diverse purposes, including river and LW management, as well as flood risk prevention.

Figure 1 .
Figure 1.Sketch of the large wood (LW) drift mechanism in a constant curvature channel due to the effect of the helicoidal secondary flow in the fully developed region.The angle δ represents the deviation of the shear stresses (τ s ) at the surface with respect to the mean flow velocity field.

Figure 2 .
Figure 2. Scheme of the implemented methodology used for simulating the effects of secondary currents in river bends in Iber-Wood.The effects on the hydrodynamics affect the large wood (LW) as the model is fully coupled (see the text for details).

Figure 3 .
Figure 3. Sketch of the considered channel bend with specifications of point gauge locations and indications about 3D flow field measurements performed by Zaid (2018), and locations of the dowel inlet points (P 1 -P 3 ) considered in the study byInnocenti et al. (2022).

Figure 4 .
Figure 4. Secondary current intensity variation along the bend (a), and the deviation angle at the free surface (b) for method M1.Negative deviation means clockwise rotation of vector direction.

Figure 5 .
Figure 5. Observed and simulated water depth (a) and surface velocity magnitude (b) for the model without any correction and the three implemented methods (M1-M3).

Figure 6 .
Figure 6.Observed and simulated flow streamlines at the water surface along the channel bend.Three flow streamlines are considered: s1, s2, and s3, respectively (a).The similarity expressed in terms of C s between observation and simulations was determined for the entire length of the bend (b), the first part of the bend (c), and the second part of the bend (d).

Figure 7 .
Figure 7. Similarity between observed and simulated wood trajectories considering different dowel inlet points and dowel initial orientations.(a, b) The results for the experimental configuration in which the dowel was inserted from P 2 with the initial orientation parallel; (c) the coefficient of similarity C s for the six considered experiments.All trajectories are related to the dowel center of mass.

Figure 8 .
Figure 8.Comparison between observed and simulated dowels travel time considering the three inlet points and the two initial dowel orientations: parallel (a) and perpendicular (b) to the channel axis.P 1 -P 3 indicate the dowels' inlet points considered in the study by Innocenti et al. (2022).

Figure 9 .
Figure 9. Sensitivity analysis of the proposed methodology (M1) to the correction term ∆c (a, b) and to the relaxation length λ (c, d).Insets report the comparison between the observed trajectory and a series of simulated trajectories (i.e., modeled trajectories for ∆c = 0.1 and ∆c = 0.8, λ calculated as in Equation 17, λ gen = λ rel = 1 m, and λ gen = λ rel = 10 m); all trajectories are related to the dowel center of mass.
Note.The reported values represent the global value of RMSE, RSE, and NSE, for water depth and water surface velocity magnitude in simulating the flow at the channel bend.The minimum and maximum measured values are reported in parentheses.Performance of the Numerical Model to Simulate the Leichtweiß-Institute Experiment: Iber Without Any Correction and Using M1-M3