Method to evaluate large-wood behavior in terms of the convection equation associated with sediment erosion and deposition

Recent flood hazards occurring in mountainous areas are often characterized by numerous amounts of sediment and large wood supplied from upstream, which often exacerbate flood disasters in downstream areas. This paper proposes a method for describing large-wood behavior in terms of the convection and storage equations, together with the governing equations for describing flood flows and channel changes associated with active sediment erosion and deposition. The proposed method is tested for its validity by simulating the phenomena occurring in an open channel with an erodible bed and flood flow with numerous amounts of sediment and large wood in the Akatani River flood disaster. As a result of calculations reproducing the open channel experiment, the applicability of the method is indicated as the percentage of wood pieces captured in the sediment deposition areas in the channel is within the range of the experimental results. The results of 2-D flood flow calculations with sediment and large wood in the Akatani River flood disaster suggested that large-wood deposition is reproduced where bed deformation is well reproduced. Overall, since the proposed method makes it possible to simulate the behavior of various amounts of large wood, it can be applied to the management of hazards in mountainous rivers such as the Akatani River.

1 Introduction

On 5 July 2017, extremely heavy rains hit northern Kyusyu, Japan, causing landslides and debris flows at numerous locations in mountainous areas. In Japan, flood disasters with a large amount of sediment deposition are becoming increasingly apparent in mountainous areas or steep slope areas. This is caused by sediment and large wood pieces transported by landslides and debris flows associated with heavy rainfall and deposited at small valley outlets, transported by flood waters, and inundated with flood waters where the sediment transport capacity is rapidly reduced.

The Akatani River, with a basin area of approximately 20 km², was particularly severely damaged in this disaster (Chen et al., 2018; Harada and Egashira, 2018). The photo on the left side of Fig. 1, taken immediately after the disaster, shows that the sediment and large wood carried by the debris flow accumulated at the valley outlet. According to the results of the laser profiler survey, approximately 2.9 million m³ of sediment was produced by landslides and debris flows in the Akatani River basin. Multiplying the area of the landslide and debris flow area by 549 m³ ha⁻¹ of timber volume per unit area, it is estimated that approximately 39 000 m³ of large wood was produced (MLIT, 2017). This means that, assuming the standing tree density to be one tree in 2 m², approximately 19 500 pieces of large wood were produced and carried to their deposited regions during the event in the Akatani River basin. Such sediment and large wood produced by landslides and debris flows are deposited in areas where the slope is about 4 to 6^∘, often forming a debris flow fan at the exit of valleys, as shown in the photo on the left side of Fig. 1. When the deposited sediment and large wood are eroded by the water flow, some is transported with the flood flow through the channel network, and once such large wood is trapped on bridges and other structures, it disrupts the flood flow. The photo on the left side of Fig. 1, taken immediately after the disaster, shows that considerable amounts of large wood accumulated at around the bridge. Recently, these types of flood hazards have been reported to occur in many steep, forested areas (Lucía et al., 2015, 2018; Steeb et al., 2017; Comiti et al., 2016; Harada et al., 2019). As reported in those studies, large wood pieces often become a major contributing factor that exacerbates flood disasters. Hence, it is highly important to develop a method to evaluate their behavior in a flood flow.

Figure 1 Large-wood deposition at the outlet of a valley bottom (left) and large-wood deposition around the bridges (right) from the Akatani River flood disaster in 2017. The left photo was taken by the authors, and the right photo was taken by the Geospatial Information Authority of Japan.

Research on large wood in rivers has been conducted, and some studies have proposed simulation models for their behavior (Swanson et al., 2021; Ruiz-Villanueva et al., 2019). Nakagawa et al. (1994) proposed a numerical simulation model to compute the behavior of individual wood pieces in a two-dimensional flow field by calculating the transportation of wood using a combination of translocation and rotation. This method was applied and verified in several laboratory experiments (Shrestha et al., 2009, 2012). Gotoh (2001) proposed a method to track the motion of driftwood based on the Lagrange particle method (MPS method) by treating wood pieces as rigid bodies. Shimizu et al. (2006) developed a model composed of two types of analyses: the Eulerian analysis of fluid motion using depth-averaged flow analysis and the Lagrange analysis of driftwood motion using the extended distinct element method. Ruiz-Villanueva et al. (2014) proposed a method to compute the transport of individual large wood pieces in a two-dimensional flow field, including physical modeling of the wood recruitment process to the flood flow. Kimura et al. (2021) developed models to compute large-wood motion in a three-dimensional flow field.

In addition to those proposed methods for the behavior of wood pieces in a flow field, some research proposed methods to analyze large-wood production processes from a watershed. Benda and Sias (2003) proposed stochastic methods to predict long-term wood budgets in watersheds. Mazzorana et al. (2009) proposed methods to estimate large-wood production processes in mountainous streams with a GIS-based index. Mazzorana et al. (2011) proposed a model for wood entrainment and deposition processes based on empirical methods and the transportation of wooden materials in the flow field. Komori et al. (2021) proposed a model to evaluate large-wood export at a watershed scale.

As suggested in the previous descriptions, the disasters that have occurred in mountainous areas in recent years can be divided into two categories: those caused by landslides and debris flows as well as their direct inundation and those caused by floods in areas where the flood and sediment transport capacity of the succeeding flood flows is rapidly reduced. The latter is generally accompanied by large amounts of sediment and large wood pieces, as flood flows erode and transport the sediment carried by landslides and debris flows. Therefore, in order to understand the characteristics of disasters in mountainous areas, it is not only important to elucidate the erosion and transport mechanisms of sediments and large wood associated with landslides and debris flows, but it is also essential to understand the erosion and transport mechanisms of sediments and large wood carried by flood flows. This study focuses on the latter topic.

Previous Lagrangian modeling attempts, such as those proposed by Shimizu et al. (2006) and Ruiz-Villanueva et al. (2014), allow for detailed evaluation of large-wood behavior by analyzing individual wood pieces in the flow. However, it is difficult to apply such Lagrangian methods to flood disasters such as the one in the Akatani River, which produced as many as 19 500 large wood pieces with active sediment erosion and deposition. Therefore, to treat such large-wood behavior with active sediment transport, this study proposes a method to describe the behavior of large wood based on a convection equation and a storage equation considering sediment erosion and deposition. Then, to investigate the validity of the proposed method in terms of erosion and deposition rates of large wood, a series of calculations are performed to reproduce the experimental results in a straight open channel. In addition, we demonstrate the application to the 2017 flood disaster in northern Kyusyu by calculating flood flows using a 2-D depth-averaged flow model with sediment and large wood. Although the reproducibility of the field application may not be perfect due to its complexity, we present it here since our target is to use the proposed method to evaluate and predict possible sediment and large-wood hazards and mitigate them in rivers.

2 Methods

To evaluate flood flows and their flooding process influenced by channel changes and large-wood traps, we employ depth-averaged 2-D governing equations consisting of mass and momentum conservation equations. The equations are expressed in the Cartesian coordinate system as follows.

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial h}{\partial t}}+{\displaystyle \frac{\partial uh}{\partial x}}+{\displaystyle \frac{\partial vh}{\partial y}}=\mathrm{0}\end{array}$$

$$\begin{array}{}\text{(2)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial hu}{\partial t}}+{\displaystyle \frac{\partial huu}{\partial x}}+{\displaystyle \frac{\partial huv}{\partial y}}=& -gh{\displaystyle \frac{\partial \left(h+z_{\mathrm{b}}\right)}{\partial x}}-{\displaystyle \frac{{\mathit{\tau }}_x}{\mathit{\rho }}}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\left({\displaystyle \frac{\partial h{\mathit{\sigma }}_{xx}}{\partial x}}+{\displaystyle \frac{\partial h{\mathit{\tau }}_{yx}}{\partial y}}\right)\end{array}\text{(3)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial hv}{\partial t}}+{\displaystyle \frac{\partial huv}{\partial x}}+{\displaystyle \frac{\partial hvv}{\partial y}}=& -gh{\displaystyle \frac{\partial \left(h+z_{\mathrm{b}}\right)}{\partial y}}-{\displaystyle \frac{{\mathit{\tau }}_y}{\mathit{\rho }}}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}\left({\displaystyle \frac{\partial h{\mathit{\tau }}_{xy}}{\partial x}}+{\displaystyle \frac{\partial h{\mathit{\sigma }}_{yy}}{\partial y}}\right)\end{array}\end{array}$$

Here, x and y are the coordinates in the major flow direction and normal to the flow direction, respectively; t is the time; h is the flow depth; u and v are the x and y components of the depth-averaged velocity; g is the acceleration due to gravity; ρ is the mass density of water; σ_xx, σ_yy, τ_xy, and τ_yx, are the depth-averaged Reynolds stresses; z_b is the bed elevation; and τ_x and τ_y are the x and y components of the bed shear stress, respectively.

Equations (1) to (3) are transformed into a general coordinate system (Shimizu and Itakura, 1991). The equations are numerically calculated by the cubic interpolated pseudo-particle (CIP) method (e.g., Yabe et al., 1991; Jang and Shimizu, 2005).

To evaluate the processes of bed deformation caused by bedload, suspended load, and suspended sediment transport in the 2-D flow field where the sediment sorting takes place actively, the following equations are employed:

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial z_{\mathrm{b}}}{\partial t}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\mathit{\lambda }}}\sum _i\left({\displaystyle \frac{\partial q_{\mathrm{b}ix}}{\partial x}}+{\displaystyle \frac{\partial q_{\mathrm{b}iy}}{\partial y}}+E_i-D_i\right)=\mathrm{0},\text{(5)}& {\displaystyle \frac{\partial c_{i}h}{\partial t}}+{\displaystyle \frac{\partial u{c}_{i}h}{\partial x}}+{\displaystyle \frac{\partial v{c}_{i}h}{\partial y}}=E_i-D_i,\end{array}$$

where λ is the porosity of the bed sediment; q_bix and q_biy are the x and y components of the bedload transport rate for grain size d_i, respectively, E_i and D_i are the erosion and deposition rates of the suspended sediment for grain size d_i, respectively; and c_i is the depth-averaged suspended sediment concentration for grain size d_i.

The bedload transport rate is estimated using a formula developed by Egashira et al. (1997), in which the constitutive relations of a water–sediment mixture flow are applied to the bedload layer:

$$\begin{array}{}\text{(6)}& q_{\mathrm{b}*i}={\displaystyle \frac{\mathrm{4}}{\mathrm{15}}}{\displaystyle \frac{K_{\mathrm{1}}^{\mathrm{2}}{K}_{\mathrm{2}}}{\sqrt{f_d+f_f}}}{\mathit{\tau }}_{*i}^{\mathrm{5}/\mathrm{2}}{p}_i,\end{array}$$

where q_b*i is the nondimensional bedload transport rate for grain size d_i; τ_*i is the nondimensional tractive force for grain size d_i; p_i is the content ratio for grain size d_i; and the other parameters, K₁, K₂, f_d, and f_f, are specified based on Egashira et al. (1997).

$$\begin{array}{}\text{(7)}& {\displaystyle }{K}_{\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{\cos \mathit{\theta }}}{\displaystyle \frac{\mathrm{1}}{\tan {\mathit{\varphi }}_{\mathrm{s}}-\tan \mathit{\theta }}}\text{(8)}& {\displaystyle }{K}_{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{‾}}{c_{\mathrm{s}}}}}{\left[\mathrm{1}-{\displaystyle \frac{h_{\mathrm{s}}}{h}}\right]}^{\mathrm{1}/\mathrm{2}}\text{(9)}& {\displaystyle }{f}_d=k_d(\mathrm{1}-e^{\mathrm{2}}){\displaystyle \frac{\mathit{\sigma }}{\mathit{\rho }}}{\stackrel{\mathrm{‾}}{c_{\mathrm{s}}}}^{\mathrm{1}/\mathrm{3}}\text{(10)}& {\displaystyle }{f}_f=k_f(\mathrm{1}-\stackrel{\mathrm{‾}}{c_{\mathrm{s}}}{)}^{\mathrm{5}/\mathrm{3}}{\stackrel{\mathrm{‾}}{c_{\mathrm{s}}}}^{-\mathrm{2}/\mathrm{3}}\end{array}$$

Here, θ is the local slope; ϕ_s is the angle of repose; $\stackrel{\mathrm{‾}}{c_{\mathrm{s}}}$ is the sediment concentration on the bedload layer; h is the water depth; σ is the density of soil particle; ρ is the density of water; e is the restitution that determines the energy loss in binary collision; k_d= 0.0828; and k_f= 0.16. Indeed, most of the values in Eqs. (7) to (10) are regarded as constant in the flow field described by Eqs. (1) to (3); thus, assuming that ϕ_s= 34^∘, $\stackrel{\mathrm{‾}}{c_{\mathrm{s}}}=\mathrm{0.25},e=\mathrm{0.85}$, and $h_{\mathrm{s}}/h\ll \mathrm{1}$, the value of the term $\mathrm{4}/\mathrm{15}{K}_{\mathrm{1}}^{\mathrm{2}}{K}_{\mathrm{2}}/\sqrt{f_d+f_f}$ in Eq. (6) is approximately 4.35.

In Eq. (8), h_s is the thickness of the bedload layer, which is described as follows (Egashira et al., 1997).

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{h_{\mathrm{s}}}{h}}={\displaystyle \frac{\mathrm{1}}{(\frac{\mathit{\sigma }}{\mathit{\rho }}-\mathrm{1})\stackrel{\mathrm{‾}}{c_{\mathrm{s}}}}}{\displaystyle \frac{\tan \mathit{\theta }}{\tan {\mathit{\varphi }}_{\mathrm{s}}-\tan \mathit{\theta }}}\end{array}$$

The grain size distribution of bed materials is evaluated based on the concept of the bedload layer, the transition layer, and the deposition layer, which was developed by Luu et al. (2006), assuming that the mass of each material is preserved.

Erosion rates E_i of suspended sediment in Eqs. (4) and (5) are evaluated using the following equations proposed by Harada et al. (2022):

$$\begin{array}{}\text{(12)}& E_i=p_i{W}_{\mathrm{e}}\stackrel{\mathrm{‾}}{c_{\mathrm{s}}},\end{array}$$

$$\begin{array}{}\text{(13a)}& {\displaystyle \frac{W_{\mathrm{e}}}{U}}={\displaystyle \frac{K}{Ri*}},\text{(13b)}& {\displaystyle }Ri*=(\mathit{\sigma }/\mathit{\rho }-\mathrm{1})\stackrel{\mathrm{‾}}{c_{\mathrm{s}}}gh/{\sqrt{u^{\mathrm{2}}+v^{\mathrm{2}}}}^{\mathrm{2}},\end{array}$$

where W_e is the entrainment velocity at the boundary between the upper water layer and the bedload layer, Ri* is the overall Richardson number, c is the depth-averaged suspended sediment concentration, and $K=\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{3}}$ (Egashira and Ashida, 1980). Note that Eqs. (6) and (12) are used taking into account the importance of sediment sorting when evaluating the flow field with active channel changes and bed deformations.

To analyze considerable amounts of large wood in the flood flow, we assume that the wood pieces behave as neutral buoyant particles that do not affect the flow structure, as this assumption allows the introduction of the convection equation. This means that large wood pieces are transported as fine particles with the flow. In addition, the exchange of these particles for large wood pieces, such as wood deposition from the water to the riverbed, wood recruitment from the riverbed to the water, and wood entrapment on structures, should also be considered. For this purpose, the convection equation is coupled with the storage equation for large wood in the channel bed, assuming that wood erosion and deposition are proportional to sediment erosion and deposition and that wood accumulation occurs at artificial structures such as bridges. This concept is consistent with the idea that when sediment is deposited on the riverbed, water is also deposited as pore water within the void of the sediment, and when sediment is eroded, pore water is also taken into the water, and this eroded sediment and pore water contain large wood, as the wood particles are treated in terms of wood concentration.

Based on these assumptions, the behavior of large wood in the flood flow is expressed using the following equations.

$$\begin{array}{}\text{(14)}& {\displaystyle }\begin{array}{rl}\partial z_{\mathrm{b}}/\partial t>\mathrm{0}:& {\displaystyle \frac{\partial c_{\mathrm{lw}}h}{\partial t}}+{\displaystyle \frac{\partial c_{\mathrm{lw}}uh}{\partial x}}+{\displaystyle \frac{\partial c_{\mathrm{lw}}vh}{\partial y}}\\ & =-c_{*}{\displaystyle \frac{\partial z_{\mathrm{b}}}{\partial t}}{c}_{\mathrm{lw}}r\left(t,x,y\right)\\ & -v_{\mathrm{n}}{c}_{\mathrm{lw}}{p}_{\mathrm{b}}\mathit{\delta }\left(x-x_i,y-y_i\right)\end{array}\text{(15)}& {\displaystyle \frac{\partial S}{\partial t}}={\displaystyle \frac{\partial z_{\mathrm{b}}}{\partial t}}{c}_{\mathrm{lw}}r\left(t,x,y\right)+v_{\mathrm{n}}{c}_{\mathrm{lw}}{p}_{\mathrm{b}}\mathit{\delta }(x-x_{i}y-y_i)\text{(16)}& {\displaystyle }\begin{array}{rl}\partial z_{\mathrm{b}}/\partial t<\mathrm{0}:& {\displaystyle \frac{\partial c_{\mathrm{lw}}h}{\partial t}}+{\displaystyle \frac{\partial c_{\mathrm{lw}}uh}{\partial x}}+{\displaystyle \frac{\partial c_{\mathrm{lw}}vh}{\partial y}}\\ & =-c_{*}{\displaystyle \frac{\partial z_{\mathrm{b}}}{\partial t}}{\displaystyle \frac{S}{D}}r\left(txy\right)\\ & -v_{\mathrm{n}}{c}_{\mathrm{lw}}{p}_{\mathrm{b}}\mathit{\delta }(x-x_{i}y-y_i)\end{array}\text{(17)}& {\displaystyle \frac{\partial S}{\partial t}}={\displaystyle \frac{\partial z_{\mathrm{b}}}{\partial t}}{\displaystyle \frac{S}{D}}r\left(t,x,y\right)+v_{\mathrm{n}}{c}_{\mathrm{lw}}{p}_{\mathrm{b}}\mathit{\delta }(x-x_{i}y-y_i)\end{array}$$

Here, c_lw is the depth-averaged large-wood concentration, c_* is the sediment concentration of the stationary bed, v_n is the inward velocity normal to the structure area such as the bridge; D is the depth of a standing tree's roots, and S is the stored volume of large wood in a unit area of the ground or the riverbed per unit area. When the volume of a piece of wood is V, S and the amount of wood (N) in a certain area (A) are converted by $S=VN/A$.

Equations (14) and (16) are the convection equations for large-wood transport with the water flow, and Eqs. (15) and (17) are the storage equations of large wood stored on the bed. The first term on the right-hand side in Eqs. (14) to (17) represents the large-wood exchange between the water flow and the channel bed. Since we assume that the erosion and deposition of large wood take place in proportion to sediment erosion and deposition, the term in Eqs. (14) and (15) represents the wood deposition from the water to the bed, and the term in Eqs. (16) and (17) represents the large-wood recruitment from the bed to the water. Figure 2a shows the concept of these processes.

Figure 2 Concept of large-wood recruitment and deposition (a) and the relation between bed erosion, root depth, and large-wood recruitment (b).

In the case of sediment deposition, $c_{*}\partial z_{\mathrm{b}}/\partial t$ corresponds to the sediment deposition rate per unit time in unit area. We assume that large wood pieces within the height range are entrained into the riverbed; thus, the amount of $c_{*}\partial z_{\mathrm{b}}/\partial t{c}_{\mathrm{lw}}$ is stored in the riverbed. In the case of sediment erosion, as shown in Fig. 2b, when the bed erosion reaches root depth D, all wood storage S is recruited to the water. Therefore, $S{c}_{*}\partial z_{\mathrm{b}}/\partial t/D$ corresponds to the large-wood recruitment from the bed to the water per unit time.

Large-wood recruitment does not occur at depths shallower than a certain water depth, and large-wood deposition does not occur at depths deeper than a certain water depth. Function r(txy) in Eqs. (14) to (17) is introduced to describe these cases and set as shown in Fig. 3 in the present research.

Figure 3 Specification of the functional form of r(txy).

In Eqs. (14) and (16), the last term defines large-wood capture at the structures. Dirac's δ function is employed to evaluate the capturing of large wood at structures such as bridges by defining the locations of these types of structures as (($x,y)=$ (x_i, y_i)) and setting Dirac's δ function for the structures as δ= 1 and for other places as δ= 0. p_b denotes the probability that large wood is captured at structures, ranging from 0 to 1.

When large wood is present in the bed, the sediment transport rate may be reduced due to the cover or shielding effect of the large wood. Therefore, in the present study, the bedload transport rate in Eq. (6) and the erosion rate of suspended sediment in Eq. (12) are reduced at the rate of the shielding effect.

At a point where a structure, such as a bridge, exists, e.g., (($x,y)=$ (x_iy_i)), the water decreases in the cross-section where the velocities across the cell, i.e., uΔx and vΔy, are reduced by uαΔx and vαΔy, in which α is described as follows:

$$\begin{array}{}\text{(18)}& \mathit{\alpha }={\displaystyle \frac{S}{h}}/(\mathrm{1}-c_{*drf}),\end{array}$$

where c_*drf is the wood concentration of a stationary layer composed of the deposited large wood only.

3 Tests of the erosion and deposition terms of large wood

3.1 Outline

In the proposed method, we assume two major factors that affect calculation results: the exchange of large-wood particles between the water and the bed as well as wood entrapment on structures. In order to investigate the validity of these factors, we compare the results of previous flume experiments and the calculations performed to reproduce the flume experiments. In the flume experiment conducted by Itoh et al. (2010), a sediment bed was prepared in the flume, and water and wood pieces were constantly supplied from the upstream end to investigate the deposition rate of the wood pieces in the flume. The original purpose of the experiment was to investigate the effect of the specific gravity of wood pieces on their deposition rate, and since the present method assumes that the erosion and deposition of large wood take place in proportion to sediment erosion and deposition, this experiment is suitable for investigating the validity of the two factors.

In some cases of the experiment where the flume width and the log length were the same, once a log was deposited in such a way as to block the cross-section of the flume, other logs were caught one by one with that first log, resulting in the log jam phenomenon. The log jam formed in the flume in run 16 is shown in Fig. 4. This log jam does not occur in the calculation without a trigger because the present method treats logs as the depth-averaged concentration of large wood. Therefore, in some cases of the calculation, obstacles that trap wood pieces, i.e., the locations where δ= 1 in Eqs. (14) to (17) with p_b=1, are set in the channel-crossing direction as a trigger so that wood is trapped there.

Figure 4 Log jam formation in run 15; the authors modified the Itoh et al. (2010) photos.

3.2 Experiments

In the flume experiment conducted by Itoh et al. (2010), a straight open channel with a length of 10 m, a width of 20 cm, and a slope of 0.045 was employed. The experimental cases and the results are shown in Table 1, where runs 1 to 8 were conducted without sediment; thus, these cases are not included in this study. In runs 9 to 15, uniform sediment with a grain size of 18.3 mm in d₆₀ was continuously supplied from the upstream end so that the sediment equilibrium condition was achieved during the experiment. The flow discharge in the experiment was kept constant at 1.42 (l s⁻¹), and pieces of wood with a length of 20 cm, equal to the width of the flume, and a diameter of 0.61 cm were randomly supplied from the upstream end for 60 s. Throughout the experiment, the log entrapment rate, i.e., the ratio of logs trapped in the flume to logs supplied, was measured. As shown in Table 1, the wood models were made of polyethylene with a specific gravity of 0.95, polymethyl methacrylate with a specific gravity of 1.20, or a 1:1 mixture of the two. The wood supply rate varied from case to case, and basically, the higher the supply rate, the higher the log entrapment rate due to the higher probability of log jam formation. In the cases where a log jam was not formed, the log entrapment rate ranged from 0 % to 3 %, and the logs were trapped where sediment deposition occurred along with the bed forms. In the cases where a log jam was formed, the log entrapment rate ranged up to 77.2 %. In cases where two log jams were formed, more logs tended to be trapped.

Table 1 Experiment cases and results from Itoh et al. (2010).

3.3 Numerical simulations

A depth-averaged 2-D flood flow model with sediment transport and bed deformation, iRIC-Nays2DH (Shimizu et al., 2019), which is partially modified by the authors, is employed to reproduce the flume experiments. The computational grid size is 1 cm², and the other computational conditions are the same as in the experiment conditions shown in Sect. 3.2. The upstream wood concentration c_lw is set to 0.41 %, which corresponds to the wood supply rate of one log per second (1 log s⁻¹) in the experiment. To induce bed deformations in the calculation, small perturbations with a grain size scale are randomly applied to the initial bed. Therefore, before the logs are supplied from the upstream end, the channel is sufficiently close to equilibrium with sand bars that have formed in the bed. The computational cases and the results of the log entrapment rates are shown in Table 2. In cases 4 to 7, obstacles that trap wood pieces are set, as shown in Fig. 5, to induce log jams.

Table 2 Calculation cases and results for the entrapment rates of wood in the channel.

Figure 5 Obstacle allocation concept in the computational domain.

In the computation case 1, sand bars are formed, as shown in Fig. 6. The dotted line in the top panel shows the front line of the sand bar just before the wood supply begins, and the solid line shows the front line of the sand bar after 60 s, meaning that sediment deposition occurs between these two lines due to the sand bar migration. As shown in Fig. 6, wood is deposited where sediment is deposited because the present method assumes that wood deposition occurs in response to sediment deposition. Function r(txy) in Eqs. (14) through (17) is expected to affect the wood deposition. As shown in Fig. 3, r(txy) is set to 0 when the ratio of water depth to wood diameter is greater than 2. This is because we assume that if there are branches in the wood, the wood will be deposited on the riverbed even if the ratio of water depth to wood diameter is greater than 1. To investigate the effect of this function, r(txy) is changed as shown in cases 2 and 3. According to the results of cases 2 and 3 in Table 2, the log entrapment rate is about 73 % greater in case 2 than the rate in case 1.

Figure 6 Relation between the computed elevation change (a) and wood deposition (b) in case 1. t= 0 s is the time just before the wood supply begins, and the solid line shows the front line of the sand bar after 60 s, meaning that sediment deposition occurs between these two lines due to the sand bar migration, and the wood deposition occurs where sediment deposition occurs.

In cases 4 to 7, the obstacle lines are set. For example, the log entrapment rate is 16.7 % in case 4, and most of the logs are caught at the obstacle. Since the log entrapment rate is proportional to the inward velocity normal to structure v_n, some logs are trapped at the obstacle, and the remaining logs flow downstream. In the experiment, not one log but several logs are deposited in the same place, which causes the formation of a jam. Case 5 is designed for this situation, with two rows of grids aligned across the channel as obstacles. As a result, 40 % more logs are captured in case 5 compared to case 4. Case 6 is designed to form two log jams, with two rows of obstacles placed in separate locations. As a result, about twice as much wood is captured in case 6 as in case 4. In case 7, two rows of obstacles are set in each jam, as in case 5. As a result, about 47 % more logs are captured in case 7 than in case 6.

4 Application to the Akatani River flood disaster

4.1 Target areas and hazard characteristics

The present method has been proposed to simulate the behavior of a large number of large wood pieces in a flood flow where active sediment transport and channel change occur. Therefore, we apply the method to the actual flood disaster in the Akatani River in 2017 to investigate the applicability and characteristics of the proposed method. The Akatani River basin is located on the right-bank side of the Chikugo River, where a large amount of sediment and large wood were produced in the 2017 flood disaster. The drainage area and the stream length of the Akatani River are approximately 20 km² and 8 km, respectively. According to Nagumo and Egashira (2019), 639 houses or buildings in the basin were damaged during the event. Figure 7 shows the Akatani River basin with debris flows and flood marks identified from aerial photos. This shows that numerous landslides and debris flows occurred in the mountainous areas, which increased damage to the downstream areas. Figure 8 compares aerial photos taken before and after the disaster. Although it is difficult to identify the river channel in the photo before the event because of its very narrow width, the photo after the event clearly shows sediment widely spreading over the valley bottom, indicating highly active sediment transport and deposition during the event.

Figure 7 The Akatani River basin with debris flows and flood marks identified from aerial photos (Nagumo et al., 2019, was modified by the authors). The background image is provided by the Geographical Information Authority of Japan. The debris flows and flood marks are color-coded to identify the tendency of sediment supply: (1) along the channels, (2) on the left-bank side of the Akatani River basin, (3) on the right-bank side of the Otoishi River basin, (4) on the right-bank side of the Akatani River basin, (5) in the Oyama River basin, (6) on the left-bank side of the Otoishi River basin, and (7) on the right-bank side of the Akatani River basin.

Figure 8 Aerial photos of the Akatani River before (left) and after (right) the flood event in July 2017. The background image is provided by the Geographical Information Authority of Japan.

Figure 9 shows the sediment size distribution observed immediately after the flood event. Longitudinal sediment sorting is clearly observed, exhibiting a tendency for the sediment size to become finer downstream. Moreover, although the average bed slope of the entire Akatani River channel is approximately $\mathrm{1}/\mathrm{70}$, the grain size of the deposited sediment is quite fine, which indicates that a large amount of fine sediment was supplied from the upstream area during the event.

Figure 9 Sediment sampling sites (left photo) and the sediment size distribution observed immediately after the flood event (right). The longitudinal distance corresponds to that of Fig. 4. The green line on the left shows the basin boundary. The red dotted line on the right is employed as the initial condition in the computation. The background image was taken from © Google Maps.

During the flood event, numerous amounts of large wood were supplied to the channels. Figure 10 shows the distribution of deposited wood length, which is identified from aerial photos taken just after the event. Three areas along the Akatani River channel and one valley in the basin shown in Fig. 7 are selected to investigate the length of each wood piece; the resolution of the aerial photos taken for these areas is sufficient for this purpose. According to Fig. 10, the distribution of wood length decreases from the inside of the valley to its outlet. Some of the large wood deposited in the outlet of the valley may have been transported downstream due to the flood flow and deposited along the Akatani River channel.

Figure 10 Distribution of wood length identified from aerial photos taken just after the event. The location of the valley is shown in Fig. 4. The identified wood length is shown as cumulative curves.

4.2 Upstream boundary conditions

To conduct a 2-D depth-averaged analysis under conditions such as the Akatani River disaster, it is necessary to evaluate the amount of sediment and large-wood inflow from the basin at the upstream boundary of the 2-D analysis area. In this study, we obtained the upstream boundary condition through an integrated method to simulate rainfall–runoff, landslide and debris flow, and sediment and large-wood transport in the river channel to obtain a time series of sediment and large wood discharged from the basin.

The rainfall–runoff–inundation (RRI) model, developed by Sayama et al. (2012), is prepared for the entire basin. The model deals with slopes and river channels separately. The flow on the slope grid cells is calculated with the 2-D diffusive wave model, while the channel flow is calculated with the 1-D diffusive wave model. On the slope grid cells, the slope stability analysis and debris flow computations (Yamazaki et al., 2016; Yamazaki and Egashira, 2019) are conducted. The occurrence of landslides is determined using the balanced equation of a force's action on an infinite slope. When landslide occurrence is detected, the surface soil in the cell is transported from the point of origin to the location where the deposition occurs using the equation of a mass system. Along with the sediment transport due to the landslide and debris flow, the standing wood there is also recruited to the debris flow and transported following the one-dimensional notation of Eqs. (14) to (17).

When the debris flow reaches the river channel grid cells, sediment and large wood are treated as sediment and large-wood supply to the river channel. In the river channel grid cells, sediment and large-wood transport is evaluated with the methods proposed by Egashira and Matsuki (2000), in which a section that includes the upstream confluence and excludes the downstream confluence point is designated as the unit channel, and the sediment and large-wood runoff for the entire basin is predicted by allocating the unit channels in series and parallel. As for the large-wood transport in the channel network, the behavior of the large wood in the unit channel follows the one-dimensional forms of Eqs. (14)–(17).

These models are applied to the Akatani River event in 2017 to estimate the time series of water, sediment, and large wood discharged from the basin. Japan Meteorological Agency analytical rainfall data are given as the rainfall data for the model. Since there is no hydrological record in the Akatani River basin, the model parameters were calibrated using data from the Kagetsu River basin, which is located east of the study basin, and the parameters are applied to estimate the flow discharge in the Akatani River basin. As a result, we estimate the peak discharge as approximately 340 m³ s⁻¹ at the 3.5 km point; the location corresponds to the upstream boundary of the 2-D flood flow computation, which is close to those of Shakti et al. (2018) (400 m³ s⁻¹) and ministry reports (MLIT, 2017) (400 m³ s⁻¹).

Parameters employed for the rainfall–runoff and landslide calculations are shown in Table 3. As for the landslide and debris flow, model parameters are validated so that sediment transport locations and total areas are close to those in Fig. 7 and ministry reports (MLIT, 2017). The initial conditions of sediment size distribution in the river channel are shown by the red dotted line in Fig. 9. As for the large-wood runoff computation from the basin, referring to our surveys and Kubota (2019), the density of standing trees is set as 0.06 (m³ m⁻²), assuming that the average diameter of a standing tree is 15 cm, the length 11.2 m, and the density per standing tree 2 m².

Table 3 Parameters employed for the rainfall–runoff and landslide computations.

The upstream boundary conditions obtained using this method are shown in Fig. 11. The figure shows the temporal variation of the basin-scale computational results for flood water, suspended sediment, and large-wood discharge at the 3.5 km point; the location corresponds to the upstream boundary of the 2-D flood flow computation. According to the figure, suspended sediment and large-wood discharge are concentrated before the flow discharge peak comes. Note that although it is difficult to validate these boundary conditions due to the complexity of the field and the lack of data, the results appear reasonable based on the information we have and are sufficient for the present purpose, which is to investigate the applicability and characteristics of the proposed method.

Figure 11 Computed results, i.e., upstream boundary conditions for 2-D computation, flood water, suspended sediment, and large-wood discharge at the 3.5 km point.

4.3 Computational conditions for the 2-D flood flow with sediment and large-wood behavior

The computation area is approximately 3.5 km long, as shown in Fig. 7. The average bed slope of the computational domain within the 3.5 km point is approximately $\mathrm{1}/\mathrm{120}$. For the computation, iRIC-Nays2DH (Shimizu et al., 2019), which was partially modified by the authors, is employed. As the initial topography, DEM data measured by an aerial laser survey before the flooding are used. The roughness coefficient is set as 0.03 for the entire computation domain, which was determined so that the flood marks would generally match the computation results. The initial sediment size distribution, indicated by the red dotted line in Fig. 9, is given within the 3.5 km reach. The grid size is 5 m by 5 m. No large-wood deposition is set; i.e., S is set to 0 in Eqs. (14)–(17) in the entire calculation domain.

Seven bridges inside the domain are set as obstacles, and δ=1 in Eqs. (14)–(17) at these locations. The large-wood capture rate p_b at bridges can take values between 0 and 1, but in this study, p_b is uniformly set to 1. In this computation, when large-wood accumulation takes place on bridges, the cross-sectional area of the flow in the grid is reduced, which in turn affects the flow conditions around the bridges. Calculations are performed for the three cases shown in Table 4 to compare the differences in results depending on the presence of sediment and large wood. Case 1 is the flow computation only without sediment and large wood, case 2 is the flow with sediment but without large wood, and case 3 is the flow with sediment and large wood.

Table 4 Calculation conditions for the 2-D flood flow with sediment and large wood.

4.4 Computation results

Figure 12 compares cases 1, 2, and 3 and shows a water depth contour map at the peak discharge time, which shows the area between 1.2 and 2.5 km from the downstream end of the computational domain. For example, the area circled by the solid white line in Fig. 9 shows that a wider area is inundated in case 3 than in case 1, which is closer to the actual inundated area.

Figure 12 Comparison of water depth at peak discharge between case 1 (left), case 2 (middle), and case 3 (right). Case 1 is the flow only, case 2 is flow with sediment but without large wood, and case 3 is flow with sediment and large wood. The white dotted line indicates the inundated area as deciphered from an aerial photo. The background image was taken from © Google Maps.

Figure 13 compares the difference between the ground elevation measured by an aerial laser survey before and after the flooding (left) and the difference between the beginning and end of the calculation for the ground elevation in case 3 (right). Although the calculation results show a little excessive sediment deposition upstream of the 1.5 km point, the two trends are generally consistent in that more than 2 m of sediment deposition takes place in the river channel and that sediment is deposited several tens of centimeters to 1 m thick in the areas where inundation occurred.

Figure 13 Comparison of elevation changes before and after the flooding measured by an aerial laser survey (left) and elevation change at the end of case 3 calculation (right). The background image was taken from © Google Maps.

Figure 14 compares the observed number of large wood pieces deposited in a 25 m² area with the number of computed pieces deposited at the end of the calculation. The number of observed large wood pieces deposited within a 25 m² area is determined from aerial photographs taken immediately after the flood event. The computed results, i.e., S in Eqs. (14)–(17), are converted to the number of large wood pieces by assuming that the diameter and length of a piece of wood are 20 cm and 7 m, respectively, referring to Fig. 10 and Kubota (2019). It should be noted that some of the deposited large wood pieces are not distinguished because they are buried.

Figure 14 Comparison of the large-wood deposition between that observed from the aerial photos (left) and computed results (right). The computed results are converted to the number of pieces of large wood by assuming that the diameter and length of a piece of wood are 20 cm and 7 m, respectively. The background image was taken from © Google Maps.

Figure 15 shows the water level in the river channel and the riverbed height at the peak flow in three different cases and compares the results with the trace water level. In case 3, the accumulation of large wood near the bridge obstructs the river channel flow; thus, the water level rises markedly upstream of the bridge. Comparing the trace water level and the calculated water level around the 1.5 km point, the water levels in cases 1 and 2 are about 1 m lower than the trace water level, and the water level in case 3 is partly due to the large-wood capture rate p_b at the bridge being uniformly set to 1, but at least the water level is evaluated lower in cases 1 and 2, where large wood is not considered. In case 3, the bed shear stress in the river channel is reduced upstream of the bridge, thus causing significant sediment deposition here.

Figure 15 Comparison of the results of each case with the water level mark and riverbed elevation in the longitudinal direction of the river channel during peak flow. The bridge location shows the elevation of the bridge height.

In terms of sediment deposition in the river channel, though the post-flood bed elevation is shown in Fig. 15, it should be noted that these are not direct comparisons, as the lines in cases 1 to 3 show the calculated bed elevation at the moment of the peak flow. The figure shows that in case 3, a large amount of sediment has already been deposited in the channel prior to the peak flow, which significantly reduces the channel capacity prior to the peak flow. Such sediment deposition is also seen in the observed post-flood elevation. In case 2, the sediment deposition in the river channel is also significant; however, the amount of sediment deposition is not as large as in case 3 because the deposition of large wood at the bridge and the associated flow obstructions are not calculated. Due to these effects, the flood inundation expands over the valley bottom, as is especially noticeable near the bridge in case 3 in Fig. 12.

Figure 16 shows the contours of the flow velocity in the vicinity of the bridge (1.2 to 1.5 km). Figure 16 compares case 2, in which large wood is not computed, with case 3, in which large wood is computed. The flow is obstructed in the bridge due to the large-wood accumulation at the bridge, causing the flow to divert around it. This results in a larger area of inundation in case 3 and a larger area subject to higher velocity fields.

Figure 16 Difference in flow pattern between case 2 (left) and case 3 (right) around the bridge at peak discharge. The background image was taken from © Google Maps.

5 Discussion

5.1 Assessments of key model assumptions

The present paper proposes a new method capable of simulating the behavior of large numbers of large wood pieces in the flow field based on the convection and storage equations. Figure 17 compares the experiment and computation results in terms of the number of log jams and the log entrapment rate. Most of the calculation results range within the experiment results, which show the validity of the present method in terms of large-wood erosion and deposition. In cases without log jams, although the log entrapment rate ranges from 0 % to 3 % within a 10 m long channel, if the channel distance is long, such as a river channel, a lot of wood will be deposited along with the flow associated with sediment deposition. In terms of wood deposition, function r(txy) is a major factor that influences the results, taking into account the diameter of the wood and also the presence of branches; in cases where the function is in case 2, the log entrapment rate is 70 % higher than in case 1. In cases with log jams, the entrapment rate of large wood, p_b, should be further investigated.

Figure 17 Comparison of the flume experiment and calculation in terms of the number of log jams and the log entrapment rate. The original data are shown in Tables 1 and 2.

In addition, in terms of the model application to the field, the accuracy of the bed deformation has a significant effect on the results. For example, Fig. 14 compares the measured and calculated results for large-wood deposition. In area (b), where large wood tends to accumulate near the bridge, the observed and calculated results are in partial agreement. On the other hand, in area (a), the observed results show that large wood is deposited far from the original river channel, while the calculated results show that large wood is deposited close to the original river channel, i.e., the right side of the white dotted rectangle. In this area, the flow that is separated from the main flow becomes an eddy and deposits suspended sediment and large wood in some places far from the original river channel, while the phenomena are not well reproduced in the computation due to the grid-scale problem. This means that in order to accurately evaluate large-wood deposition within a 2-D flow scale, a fine mesh must be used to evaluate sediment deposition.

5.2 Comparisons to previous modeling attempts

Previous Lagrangian modeling attempts, such as those proposed by Shimizu et al. (2006) and Ruiz-Villanueva et al. (2014), can accurately analyze the behavior of large wood by tracking individual pieces of wood. On the other hand, in order to evaluate the behavior of large wood in disasters such as the Akatani River, where as many as 19 500 pieces of large wood were produced, it is difficult to use such existing Lagrangian methods, while the present method can overcome such difficulties. The present method also allows for basin-scale analysis of the production, transport, and deposition processes of large wood, which is useful even when the basin covers a large area. In a 2-D flow model, large-wood behavior is easily treated even in a field with significant riverbed deformations.

On the other hand, as already described, the disadvantage of this approach is that the conversion from large-wood concentration to actual large wood pieces and vice versa, such as wood deposition from the water to bed, wood recruitment from the bed to water, and large-wood entrapment on structures, is not necessarily adequate at this stage. These points can be improved in the future by performing hydraulic experiments or, for example, by using the method of Kimura et al. (2021), which allows accurate calculations of the behavior of large wood in a flow field with wood branches. Another disadvantage, since large-wood erosion and deposition is proportional to the sediment erosion and deposition, is that the precision of the bed deformation calculation directly affects the results of large-wood erosion and deposition. In this sense, in order to obtain a reasonable result, it is necessary to use a sufficiently fine mesh when computing a 2-D flow model that can reproduce, for example, an eddy separated from a main flow.

6 Conclusions

The present paper has proposed a method to evaluate the behavior of large wood in the flow field based on the convection equation and the storage equation with active sediment transportation and channel bed deformation, which characterizes recent flood disasters in mountainous and hilly regions, such as the flood disaster in the Akatani River in 2017. Within the method, to investigate the validity of the erosion and deposition terms of large wood, flume experiment results are compared with the calculation results. The experimental results show that large wood is trapped in sediment deposition areas, which are also reproduced by numerical simulations. In terms of the wood deposition, the function $r\left(t,x,y\right)$ is a major factor that influences the results; varying it within realistic values resulted in a maximum difference of about 70 % in the rate of large-wood deposition, and these results are still in the range of the experimental results. In cases with log jams forming, though the large-wood capture rate p_b should be further investigated, the large-wood capture rates are within the range of the experimental results, indicating the validity of the proposed method. The results of 2-D flood flow calculations with sediment and large wood in the Akatani River flood disaster suggested that large-wood deposition is reproduced where bed deformation is well reproduced. However, there are areas in the computational results where the flow pattern and subsequent bed deformation are not properly reproduced, and the large-wood deposition is not well reproduced.

Since the proposed method makes it possible to simulate the behavior of a large number of large wood pieces, it can be applied to the management of hazards, such as flooding in the Akatani River. The computed results are useful to assess the effectiveness of countermeasures, develop hazard maps, and create evacuation plans. In addition, the effect of countermeasures such as large-wood-capturing structures can be evaluated through simulations using the proposed method, which provides practical information to control hazards more efficiently and effectively.