Numerical Study on the Shear Stress Characteristics of Open-Channel Flow over Rough Beds

: Bed shear stress is an important measure of benthic habitats since it is related to many ecological processes. In this study, we focused on the ﬂuctuating characteristics of shear stress in rough-bed open-channel ﬂows. The roughness element method was adopted to mimic natural rough beds and the Improved Delayed Detached Eddy Simulation (IDDES) model was used to obtain comprehensive information about shear stress near the rough bed. Three arrangement patterns of the roughness elements were simulated to compare their effects on ﬂow structure and shear stresses. The arrangements of the roughness elements altered the Reynold stress and turbulent kinetic energy characteristics, due to the variance of blockage in lateral directions that led to ﬂow detachment and changes in the ﬂow directions. Quadrant analysis revealed the spatial variations of the instantaneous shear stress burst events at different locations in the wake. By using spectrum analysis, the accumulation of shear-stress energy from small to large vortex scales was estimated, which revealed that the instantaneous effect of the shear stress was signiﬁcantly stronger than the effect of the time-averaged shear stress, especially on small scales. The results of this study suggest the signiﬁcance of the ﬂuctuation part of shear stress in further studies on ecological processes.


Introduction
The construction and operation of dams and hydropower plants dramatically change river hydrology and hydrodynamics, involving a series of hydrodynamic variables such as flow velocity, depth, and shear stress [1,2]. Among these variables, shear stress is an important measure of benthic habitats since it plays a critical role in many ecological processes including macroinvertebrate drifting [3] and bio-particle transportation [4]. An increase in shear stress destabilizes benthic macroinvertebrate habitats by moving sediments or directly forcing benthic macroinvertebrates to leave their habitats, which both lead to spatial variations in the benthic macroinvertebrate community [5][6][7]. Additionally, shear stress is widely regarded as the key indicator of bed material movements, which makes it significant in studies of river morphology.
In natural rivers, shear stress is a fluctuating signal and its related ecological processes may not only depend on its time-averaged value (usually regarded as the Reynolds stress) but also on its fluctuations. Some studies have highlighted the significance of the instantaneous fluctuation part of shear stress in turbulence [8][9][10]. However, the fluctuation part of shear stress has been seldom studied, compared to the time-averaged shear stress studied widely [5,11,12]. Known as 'burst events', shear stress fluctuations can be divided into four event types by quadrant analysis, among which ejections and sweeps are dominant. The fluctuation part of shear stress generates an additional force on the particles and thus reduces the shear stress threshold of dislodgement [9,10,13]. Blanckaert et al. revealed that the temporal fluctuation of flow velocity results in an increase in drag force by about 100% [9]. Turbulence fluctuation is related to the development of vortex structures, and the magnitude of fluctuation is significantly affected by vortex scales. It is indicated that the magnitude of shear stress may vary with the vortex scale, although the spectral analysis of turbulence is rarely adapted to shear stress.
The direct measurement of bed shear stress is a difficult task, and its spatial and temporal characteristics are usually analyzed based on near-bed hydrodynamic conditions. Flow field data can be acquired by direct measurements, such as Particle Image Velocimetry (PIV), Acoustic Doppler Velocity Profiler (ADVP), and intrusive Electromagnetic Current Meters (ECMs) [8,9,14,15]. However, these high-resolution measurements that satisfy the analysis of local shear stresses have a high cost and the results are also affected by the wall boundary or measurement instrument, and by now near-bed flow structures are still poorly investigated [16]. An alternative method is numerical simulation, by which the disadvantages in measurements could be overcome, and the comprehensive flow field data are available. Currently, in the ecohydraulic field, the most commonly used hydrodynamic models (e.g., PHABSIM, RHYHABISIM, and River-2D) only involve timeand space-averaged flow variables [17][18][19][20]. Morris et al. tried a three-dimensional Large Eddy Simulation (LES) in their study on Glossosoma biomass distribution, based on a high-resolution riverbed morphology measurement [21]. This provides insight into the application of three-dimensional high-resolution flow simulations in the ecohydraulic field, but this costly case study cannot be used widely. Thus, an adequately simplified method of rough channel beds is needed to balance computation costs and accuracy.
In numerical simulations, near-bed flows are commonly modeled by three methods, namely wall functions, wall with roughness elements, and wall applying measured rough bed morphology [21][22][23]. Wall functions over-simplify the rough bed and are unable to reflect the effect of large-size sediments on the local flow field. Additionally, as mentioned above, rough bed morphology measurements are costly. Roughness elements can be an appropriate tool to capture the flow characteristics of the roughness of a natural riverbed [16]. Among the many types of roughness elements, two-dimensional bars and three-dimensional cubes are most frequently used. Previous numerical studies have revealed the influence of the roughness element dimension, spacing ratios between the bar interval p and roughness element height d, cross-section shapes, and arrangement patterns upon turbulence characteristics such as turbulence intensities and drag coefficients [22][23][24][25][26]. In comparison, numerical studies and statistical analyses on shear stress characteristics, especially the fluctuation part of shear stress, are seldom found.
The main purpose of this study is to evaluate the shear stress characteristics in turbulent flow, especially the instantaneous characteristics, based on numerical results by using the roughness element method to mimic riverbeds. We first introduce the turbulent model applied in this study and test its validity, and then describe the roughness element arrangements of the conducted simulations. Following are the results of the numerical simulations analyzed by quadrant analysis and spectrum analysis. Discussions about the shear stress fluctuation events' spatial variation and the accumulated effect of shear stress at different scales are also included in this section. Finally, the concluding remarks are presented.

Improved Delayed Detached Eddy Simulation Model
This study used the Improved Delayed Detached Eddy Simulation (IDDES) model to investigate the turbulent characteristics of near-bed flow, considering the balance between calculation cost and accuracy. The IDDES model is a Reynolds-averaged Navier-Stokes (RANS)/Large Eddy Simulation (LES) hybrid model based on the Detached Eddy Simulation (DES) model, whose feasibility in simulating turbulent fluctuations has been proved by previous research [27]. The model invokes either the RANS model or LES model in  different flow regions according to the relative size of the turbulent length scale and grid  scale. The regions in which the turbulent length scale exceeds the grid scale are solved  by the LES model, and the regions near the solid boundaries and the locations in which  the turbulent length scale is smaller than the maximum grid scale are solved by (1) and (2): in which ρ is the fluid density, k is the turbulent kinetic energy, ω is the specific turbulent dissipation rate, → U is the flow velocity vector, µ is the fluid viscosity, µ t is the turbulent eddy viscosity, P k is the rate of production of k, and F 1 is the blending function for α and β, α = α 1 ·F 1 + α 2 ·(1 − F 1 ), in which α 1 = 5/9 and α 2 = 0.44, and β = β 1 ·F 1 + β 2 ·(1 − F 1 ), in which β 1 = 0.075 and β 2 = 0.828. l IDDES is the IDDES length scale. σ k = 2. σ ω2 = 0.856 [28].
The original DES approach will incorrectly invoke the LES model when the grid is locally refined in multiple directions, especially in boundary regions. As a result, the eddy viscosity may be substantially underestimated from the RANS region to the LES region, triggering Grid Induced Separation. In the IDDES model, intricate blending and shielding functions are applied in characteristic length scales to protect the wall boundary layer from the incorrect invocation of the RANS model [28]. Equation (3) shows the length scale used in this model: in which l RANS , the RANS turbulence length scale, is calculated from the turbulence kinetic energy k and the specific rate of dissipation ω for the k − ω model, i.e., l LES is the LES filter width, defined as: with C DES = 0.65, C w = 0.15 , and ∆ max = max(∆x, ∆y, ∆z).f d is the empirical blending function that controls the selection of the characteristic length scale (and consequently the model applied in the grid) and f e is the elevating function that prevents an excessive reduction in the eddy viscosity in the vicinity of the RANS/LES interface. The details off d and f e can be found in Gritskevich et al., 2012. [28]. When the near-wall mesh is fine enough, the region using RANS would be very limited; namely, the majority of the domain would be calculated by LES. In this study, we adopted the SST-IDDES model and the RANS region thickness was about 0.3 mm, which was thin enough to be neglected compared with the roughness height and the total flow depth.
The simulations were conducted by using Ansys Fluent 14.0, where the governing equations are numerically solved based on the finite volume method and PISO scheme.

Model Validation
The model was validated by predicting a pressure-driven flow in a channel with roughness elements and comparing the prediction against the experimental results found by Krogstad, et al. [29].
The plane channel in Krogstad's experiment was 5 m long, 1.35 m wide, and 0.1 m high, with the bottom and top walls equipped with transverse square rods. The height of the rods was d = 1.7 mm and the distance between the front face of two adjacent rods was p = 8d. The cross-section average velocity was 0.1208 m/s and the Reynolds number based on the shear velocity u * was Re * = ρu * h/µ = 600, in which h was the half-height of the channel.
The bottom half of the plane channel with 15 rods was simulated in the present study, and the computational domain was L x = 0.204 m long and L y = 0.05 m high, with the roughness elements arranged on the bottom surface. The symmetry boundary condition was specified on the top boundary, and the periodic boundary condition was applied on the upstream and downstream boundaries. The width of the domain was L z = 0.2 m and the periodic boundary condition was applied to the boundaries in the lateral direction.
The calculation domain and grid of the numerical model are shown in Figure 1. The computational domain was divided into 3 layers of different resolutions in the vertical direction in order to reduce the total grid amount while maintaining enough grid resolution in the near-wall region and a good aspect ratio in the domain. The layers were connected by the interface boundary condition. The total grid amount was 789,000. The grid size in the layer near the rough floor was 0.34 × 0.34 × 0.34 mm (∆y + ≈ 4.1) while the grid size near the upper boundary was 1.36 × 3 × 3 mm. In the vertical direction, the cell height increased with a growth rate of 1.2. x, y, and z were the longitudinal, vertical, and lateral directions. u, v, and w were, respectively, the velocity components in these three directions, and u , v , and w were the fluctuation parts of the velocities. u v represented the flow shear stress exerted by the velocity fluctuations, and its time-averaged value multiplied by the flow density ρ is the Reynolds shear stress (usually represented by −ρu v ).
The model was validated by predicting a pressure-driven flow in a channel with roughness elements and comparing the prediction against the experimental results found by Krogstad, et al. [29].
The plane channel in Krogstad's experiment was 5 m long, 1.35 m wide, and 0.1 m high, with the bottom and top walls equipped with transverse square rods. The height of the rods was d = 1.7 mm and the distance between the front face of two adjacent rods was p = 8d. The cross-section average velocity was 0.1208 m/s and the Reynolds number based on the shear velocity * was * = * ℎ/ = 600, in which h was the half-height of the channel.
The bottom half of the plane channel with 15 rods was simulated in the present study, and the computational domain was Lx = 0.204 m long and Ly = 0.05 m high, with the roughness elements arranged on the bottom surface. The symmetry boundary condition was specified on the top boundary, and the periodic boundary condition was applied on the upstream and downstream boundaries. The width of the domain was Lz = 0.2 m and the periodic boundary condition was applied to the boundaries in the lateral direction.
The calculation domain and grid of the numerical model are shown in Figure 1. The computational domain was divided into 3 layers of different resolutions in the vertical direction in order to reduce the total grid amount while maintaining enough grid resolution in the near-wall region and a good aspect ratio in the domain. The layers were connected by the interface boundary condition. The total grid amount was 789,000. The grid size in the layer near the rough floor was 0.34 × 0.34 × 0.34 mm (Δy + ≈ 4.1) while the grid size near the upper boundary was 1.36 × 3 × 3 mm. In the vertical direction, the cell height increased with a growth rate of 1.2. x, y, and z were the longitudinal, vertical, and lateral directions. u, v, and w were, respectively, the velocity components in these three directions, and u', v', and w' were the fluctuation parts of the velocities. u'v' represented the flow shear stress exerted by the velocity fluctuations, and its time-averaged value multiplied by the flow density ρ is the Reynolds shear stress (usually represented by − ). Previous studies have concluded from the power spectrum that the major part of the water flow frequency ranges from 0-100 Hz. Considering the Fourier transform used in this research to analyze the turbulence characteristics, the sampling frequency should be no less than 200 Hz [30,31]. In this study, the time step and sampling interval were set to be 0.005 s, and the total sample size N = 6000 [30]. Previous studies have concluded from the power spectrum that the major part of the water flow frequency ranges from 0-100 Hz. Considering the Fourier transform used in this research to analyze the turbulence characteristics, the sampling frequency should be no less than 200 Hz [30,31]. In this study, the time step and sampling interval were set to be 0.005 s, and the total sample size N = 6000 [30].

Validation Results
The normalized velocity, turbulence intensities, and Reynolds stress of the numerical results were compared with Krogstad et al.'s measurement for model validation [29]. Their profiles were averaged horizontally over the computational domain. Figure 2 shows the velocity profiles which were normalized by the shear velocity u * .
Based on: in open-channel turbulence, in which H represents the flow depth (and equals to h here), as −u v ν ∂u ∂y in the outer region, u * can be calculated from the intercept of the linear part of −u v ∼ y at y = 0. The non-dimensional quantities are marked by the '+' superscript in this paper. The velocity's relative error of the IDDES model was calculated on each point on the profile, and their average was 5.0%. The normalized turbulence intensities and Reynolds shear stress of the IDDES model and Krogstad, et al.'s measurement are shown in Figure 3, where the turbulence intensities are represented by the root mean square error of the flow velocities and are marked by the 'rmse' subscript. The average relative errors of the turbulence intensity profiles in the x, y, and z directions were 9.5%, Water 2022, 14, 1752 5 of 18 9.0%, and 4.3%, respectively. As for the Reynolds stress profile, the absolute error was calculated, because the normalized Reynolds stress value approached zero on the top boundary. The absolute average error of the normalized Reynolds stress was 0.151, which was especially higher in the near-wall region where the numerical results exceeded the measurements. This was also observed in Krogstad, et al.'s direct numerical simulation (DNS) results, as the intrusive hot-wire method applied in the research caused errors in the vertical velocity measurements in the near-wall region, which tended to underestimate the measured Reynolds stress peak compared with the DNS. In other words, the real peak of the normalized shear stress should have been higher than the measured results in Krogstad, et al. [29]. For y/h > 0.2, the average error of the numerical model was reduced to 0.028, indicating that the agreement of the IDDES model was satisfactory. The normalized velocity, turbulence intensities, and Reynolds stress of the numerical results were compared with Krogstad et al.'s measurement for model validation [29]. Their profiles were averaged horizontally over the computational domain. Figure 2 shows the velocity profiles which were normalized by the shear velocity * .  Figure 3, where the turbulence intensities are represented by the root mean square error of the flow velocities and are marked by the 'rmse' subscript. The average relative errors of the turbulence intensity profiles in the x, y, and z directions were 9.5%, 9.0%, and 4.3%, respectively. As for the Reynolds stress profile, the absolute error was calculated, because the normalized Reynolds stress value approached zero on the top boundary. The absolute average error of the normalized Reynolds stress was 0.151, which was especially higher in the near-wall region where the numerical results exceeded the measurements. This was also observed in Krogstad, et al.'s direct numerical simulation (DNS) results, as the intrusive hot-wire method applied in the research caused errors in the vertical velocity measurements in the near-wall region, which tended to underestimate the measured Reynolds stress peak compared with the DNS. In other words, the real peak of the normalized shear stress should have been higher than the measured results in Krogstad, et al. [29]. For y/h > 0.2, the average error of the numerical model was reduced to 0.028, indicating that the agreement of the IDDES model was satisfactory. It can also be seen from Figure 3 that the specification of the symmetry boundary on the upper boundary will decrease the vertical turbulent intensity and increase the lateral turbulent intensity. However, the influence is limited in the region of y/h > 0.9. In the region far from the upper boundary, especially in the region near the bottom floor (where we are mostly concerned), the influence of the symmetry boundary is neglectable. Therefore, it can be assumed that, in this study, the symmetry boundary could be used to rep- For a uniform flow, the usage of periodic boundary conditions on upstream and downstream boundaries is valid for time-averaged velocities but not for velocity fluctuations. To evaluate the influence of the periodic boundary condition, spatial correlation coefficients between the inlet and the downstream points on the centerline at y = 1 mm were calculated and are shown in Figure 4. It was observed that the use of the periodic boundary condition only affected a small range near the boundaries, indicating that in the It can also be seen from Figure 3 that the specification of the symmetry boundary on the upper boundary will decrease the vertical turbulent intensity and increase the lateral Water 2022, 14, 1752 6 of 18 turbulent intensity. However, the influence is limited in the region of y/h > 0.9. In the region far from the upper boundary, especially in the region near the bottom floor (where we are mostly concerned), the influence of the symmetry boundary is neglectable. Therefore, it can be assumed that, in this study, the symmetry boundary could be used to represent the free water surface.
For a uniform flow, the usage of periodic boundary conditions on upstream and downstream boundaries is valid for time-averaged velocities but not for velocity fluctuations. To evaluate the influence of the periodic boundary condition, spatial correlation coefficients between the inlet and the downstream points on the centerline at y = 1 mm were calculated and are shown in Figure 4. It was observed that the use of the periodic boundary condition only affected a small range near the boundaries, indicating that in the near-wall region, the longitudinal correlation of the velocity fluctuation was weak, and that the periodic boundary condition had a minor effect on the near-wall flow field. This enables the application of the periodic boundary condition for the purpose of reducing the total grid amount and improving the efficiency of the simulation. For a uniform flow, the usage of periodic boundary conditions on upstream and downstream boundaries is valid for time-averaged velocities but not for velocity fluctuations. To evaluate the influence of the periodic boundary condition, spatial correlation coefficients between the inlet and the downstream points on the centerline at y = 1 mm were calculated and are shown in Figure 4. It was observed that the use of the periodic boundary condition only affected a small range near the boundaries, indicating that in the near-wall region, the longitudinal correlation of the velocity fluctuation was weak, and that the periodic boundary condition had a minor effect on the near-wall flow field. This enables the application of the periodic boundary condition for the purpose of reducing the total grid amount and improving the efficiency of the simulation.

Simulation Domain Description and Rough Wall Configuration
In Krogstad et al.'s case, the half-channel height was h = 0.05 m and the roughness height d = 1.7 mm, which is much lower than the flow depth and roughness height caused by substrates in natural rivers, respectively. Therefore, the following simulations were based on Blanckaert et al.'s experiment to approximate a turbulent flow over a natural rough bed [9].

Simulation Domain Description and Rough Wall Configuration
In Krogstad et al.'s case, the half-channel height was h = 0.05 m and the roughness height d = 1.7 mm, which is much lower than the flow depth and roughness height caused by substrates in natural rivers, respectively. Therefore, the following simulations were based on Blanckaert et al.'s experiment to approximate a turbulent flow over a natural rough bed [9].
The open channel used in Blanckaert et al.'s experiments were 8.5 m long and 0.5 m wide. We adopted one test from their experiments, in which the flow depth was 0.2 m and the average flow velocity was 0.3 m/s. The bottom sediment characteristic sizes were D 50 = 0.8 mm, D m = 2.3 mm, and D 90 = 5.7 mm [9]. In the following simulations, the computational domain was L x = 0.24 m long, L y = 0.2 m high, and L z = 0.24 m wide. The boundary conditions were the same as those in the validation case. The three-layer layout of the mesh was also adopted, with the total grid amount ranging from 3,474,600 to 6,568,800 depending on specific cases. The grid size near the rough floor was about 0.5 mm × 0.5 mm × 0.5 mm, while the grid size near the water surface was 2 mm × 4 mm × 2 mm. The time step ∆t and sample size were identical to the validation case. Sampling started after a simulation duration of 120 s when turbulence was fully developed in the domain.
The roughness elements were arranged on the bottom surface to mimic rough beds. Three scenarios with different roughness elements were investigated in the present study. They varied in terms of dimension, shape, spacing, and arrangements: Case 1: two-dimensional bars with square cross-section shapes were perpendicularly arranged to the flow direction ( Figure 5a). a simulation duration of 120 s when turbulence was fully developed in the domain.
The roughness elements were arranged on the bottom surface to mimic rough beds. Three scenarios with different roughness elements were investigated in the present study. They varied in terms of dimension, shape, spacing, and arrangements: Case 1: two-dimensional bars with square cross-section shapes were perpendicularly arranged to the flow direction (Figure 5a). Case 2: three-dimensional roughness elements of cubes were arranged in chessboard formation (Figure 5b).
Case 3: three-dimensional cubes were staggeringly arranged ( Figure 5c). The spacing between the roughness elements in the longitudinal and lateral directions were larger than that of the chessboard formation, letting the fluid flow by their sides. The main difference between Case 2 and 3 was that in Case 2, no interstice that flow could go through existed between the cube elements connected to each other at the corners, while in Case 3, the spaces between the cubes were wide enough to let the flow go around the cubes.
We used these models to mimic the natural riverbeds where macroinvertebrates usually shelter themselves among substrates; therefore, the interstices and roughness element height were determined according to the benthic macroinvertebrate body size (1 mm to several centimeters). In the present study, the interstices were set to be 30 mm (except for the 3D-cu-s, in which the interstices were intentionally expanded). The roughness elements height was set to be d=8 mm. Note that, by definition, the distance between the roughness elements p refer to the distance of a roughness element from its successive one Case 2: three-dimensional roughness elements of cubes were arranged in chessboard formation (Figure 5b).
Case 3: three-dimensional cubes were staggeringly arranged ( Figure 5c). The spacing between the roughness elements in the longitudinal and lateral directions were larger than that of the chessboard formation, letting the fluid flow by their sides. The main difference between Case 2 and 3 was that in Case 2, no interstice that flow could go through existed between the cube elements connected to each other at the corners, while in Case 3, the spaces between the cubes were wide enough to let the flow go around the cubes.
We used these models to mimic the natural riverbeds where macroinvertebrates usually shelter themselves among substrates; therefore, the interstices and roughness element height were determined according to the benthic macroinvertebrate body size (1 mm to several centimeters). In the present study, the interstices were set to be 30 mm (except for the 3D-cu-s, in which the interstices were intentionally expanded). The roughness elements height was set to be d = 8 mm. Note that, by definition, the distance between the roughness elements p refer to the distance of a roughness element from its successive one in the x or z direction, and p varies according to the roughness element arrangements in different cases. Figure 5 labels the definition of p in each case. Specific magnitudes of roughness element size and spacing, as well as their arrangement patterns, are listed in Table 1. These cases were named by their dimension, roughness element shape, and arrangement patterns, and are also listed in Table 1. The main characteristic hydrodynamic variables include shear velocity u * , the roughness function ∆U + , and the roughness equivalent height k s . In this study, the shear velocity u * was calculated using Equation (6). The roughness function ∆U + was calculated according to the log profile of the time-averaged longitudinal flow velocity: where the Karman constant κ = 0.41, B = 5.5. On rough surfaces, the theoretical datum y = y 0 should be considered for the application of the log profile of flow velocity. In the following part of this study, we redefine y + = 0 to be acquired at the theoretical datum y = y 0 , giving: and y + = y−y 0 ν/u * . k s is derived by: In common practice, k s /d is used for analysis.

Quadrant Analysis Methodology
Quadrant analysis was first conducted by Wallace et al. to explore the unknown information on Reynolds shear stresses from velocity fluctuations [32]. The principle of this method is to classify shear stress u v into four categories by instantaneous turbulent velocity fluctuations (u , v ): Q1 (+u , +v ), Q2 (−u , +v ), Q3 (−u , −v ), and Q4 (+u , −v ), each of which represents an event of the turbulent burst phenomenon. In this way, shear stress characteristics including magnitude, contributions of each event, and velocity fluctuation components can be investigated by quadrants [9,33].
In this study, the velocity fluctuations were locally affected by the wakes behind the different roughness elements. At different locations of the wakes, the instantaneous shear stress −u v was plotted by its u and v components on the u -v plane to visualize the shear stress characteristics, such as the magnitude of shear stress, its velocity fluctuation components, and the quadrant it is in. The magnitude of the shear stress of each point is represented by the product of its coordinate values, so points with the same absolute magnitude of shear stress are located at the same |u v | = const hyperbolas.
Previous studies suggested a spatial variation in the time-averaged shear stress on the horizontal plane near the rough floor, and that its distribution depends on the pattern of the roughness element arrangements [16]. It is necessary to classify the locations around the roughness elements and to conduct an analysis by different locations. Considering the periodic arranging patterns of the roughness elements in the present study, the locations over the roughness elements and the locations behind the roughness elements were selected since they were, respectively, under the control of the main flow and the recirculation flow in the cavities. In the quadrant analysis of the instantaneous shear stress, we chose locations above each roughness element at y = 10 mm and 9 mm behind the front edge of the roughness element (type A) and locations 18 mm behind the back edge of each roughness element in the cavities at y = 6 mm (type B). This selection is shown in Figure 6 and would be applied in all cases for further comparison. The instantaneous shear stress in the sampled 6000 steps was analyzed. lected since they were, respectively, under the control of the main flow and the recircula-tion flow in the cavities. In the quadrant analysis of the instantaneous shear stress, we chose locations above each roughness element at y = 10 mm and 9 mm behind the front edge of the roughness element (type A) and locations 18 mm behind the back edge of each roughness element in the cavities at y = 6 mm (type B). This selection is shown in Figure 6 and would be applied in all cases for further comparison. The instantaneous shear stress in the sampled 6000 steps was analyzed. In common quadrant analysis practices, longitudinal and vertical velocity fluctuations are considered as the studied flow is mainly longitudinal and lateral flow velocity could be neglected. In this study, for the type B locations, the roughness elements led to a strong three-dimensional flow in the roughness element layer, so the longitudinal and lateral velocity could be combined into a horizontal flow velocity.

Frequency Spectrum Analysis for Shear Stress
In general, spectrum analysis has been widely used to find the dominant frequency, focusing on the intensity of a vortex of a certain frequency (or certain length scale). Commonly analyzed variables in spectrum analysis include velocity, pressure, and turbulent kinetic energy [34,35].
Spectrum analysis was performed by using Fourier transform, by which a temporal signal A(t) can be transformed into a combination of trigonometric signals: In Equation (10), the constant term (the zeroth-degree harmonic term) represents the time average of the signal, and the trigonometric terms represent the signal fluctuation around the average, due to vortices of all sizes. This separates the effects of vortices of different frequencies, i.e., of different length scales.
Fourier transform is usually applied to the instantaneous velocity ( ) and the total energy of this Fourier series corresponds to the total energy of the kinetic energy of the flow. According to the Parseval theorem, this energy is the integral of the energy contained in the signal of each frequency and could be calculated from the amplitude of each harmonic term: In the Fourier series of instantaneous velocity, 2 ⁄ = , ∑ ( + ) = , therefore the integral should be + 2 ⁄ (the symbols with an overbar represent the timeaveraged values). This integral is the sum of the effects of all sizes of vortices.
In this study, we followed this idea by analogizing the fluctuating part of shear stress to the fluctuating velocity and accumulating the effects of vortices from low frequency to high frequency on shear stress. In this way, we discovered the differences between the same turbulent shear stress acting upon particles such as benthic macroinvertebrates of different sizes. Similarly, Fourier transform was applied to the instantaneous shear stress In common quadrant analysis practices, longitudinal and vertical velocity fluctuations are considered as the studied flow is mainly longitudinal and lateral flow velocity could be neglected. In this study, for the type B locations, the roughness elements led to a strong three-dimensional flow in the roughness element layer, so the longitudinal and lateral velocity could be combined into a horizontal flow velocity.

Frequency Spectrum Analysis for Shear Stress
In general, spectrum analysis has been widely used to find the dominant frequency, focusing on the intensity of a vortex of a certain frequency (or certain length scale). Commonly analyzed variables in spectrum analysis include velocity, pressure, and turbulent kinetic energy [34,35].
Spectrum analysis was performed by using Fourier transform, by which a temporal signal A(t) can be transformed into a combination of trigonometric signals: In Equation (10), the constant term (the zeroth-degree harmonic term) represents the time average of the signal, and the trigonometric terms represent the signal fluctuation around the average, due to vortices of all sizes. This separates the effects of vortices of different frequencies, i.e., of different length scales.
Fourier transform is usually applied to the instantaneous velocity u(t) and the total energy of this Fourier series corresponds to the total energy of the kinetic energy of the flow. According to the Parseval theorem, this energy is the integral of the energy contained in the signal of each frequency and could be calculated from the amplitude of each harmonic term: In the Fourier series of instantaneous velocity, a 0 /2 = u, ∑ ∞ n=1 a 2 n + b 2 n = u 2 , therefore the integral should be u 2 + u 2 /2 (the symbols with an overbar represent the timeaveraged values). This integral is the sum of the effects of all sizes of vortices.
In this study, we followed this idea by analogizing the fluctuating part of shear stress to the fluctuating velocity and accumulating the effects of vortices from low frequency to high frequency on shear stress. In this way, we discovered the differences between the same turbulent shear stress acting upon particles such as benthic macroinvertebrates of different sizes. Similarly, Fourier transform was applied to the instantaneous shear stress signal and its energy was analyzed. The accumulation of the shear stress signal energy from low frequency to high frequency could be written as a function of frequency: where f is the frequency of the shear stress signal, a 2 0 /4 = u v 2 , and ∑ f n=1 a 2 n + b 2 n = u v − u v 2 when f→∞. Compared to the frequency, the vortex length scale is a more practical indicator when relating hydrodynamic characteristics to specific objects in the flow.
It is natural to transfer the frequency value into a vortex length scale value by adopting the local time-averaged flow velocity u c as the characteristic velocity: An instantaneous shear stress signal of type B was selected for spectrum analysis to reflect the shear stress characteristics in the roughness element layer. Figure 7 displays the predicted and measured velocity profiles, normalized by the shear velocity u * (calculated by Equation (6) and listed in Table 2). These profiles have approximately equal slopes and their different intercepts separate them from each other. The roughness functions ∆U + were obtained by subtracting the intercepts from 5.5 and are listed in the third row of Table 2. Then, the roughness equivalent heights k s /d were derived by using Equation (9) and are listed in the fourth row of Table 2.

Characteristics of Time-Averaged Hydrodynamic Parameters
( ) = ( = 0) where f is the frequency of the shear stress signal, 4 ⁄ = ( ) , and ∑ ( + ) = ' ' − ' ' when f→∞. Compared to the frequency, the vortex length scale is a more practical indicator when relating hydrodynamic characteristics to specific objects in the flow. It is natural to transfer the frequency value into a vortex length scale value by adopting the local time-averaged flow velocity uc as the characteristic velocity: An instantaneous shear stress signal of type B was selected for spectrum analysis to reflect the shear stress characteristics in the roughness element layer. Figure 7 displays the predicted and measured velocity profiles, normalized by the shear velocity * (calculated by Equation (6) and listed in Table 2). These profiles have approximately equal slopes and their different intercepts separate them from each other. The roughness functions Δ were obtained by subtracting the intercepts from 5.5 and are listed in the third row of Table 2. Then, the roughness equivalent heights ks/d were derived by using Equation (9) and are listed in the fourth row of Table 2.  With the same flow rate, the roughness element arrangements had a great impact on the flow resistances. According to the values of ks/d and Δ , 2D-sq generated the largest flow resistance, followed by 3D-cu-c and 3D-cu-s, which corresponds to Volino et al.'s conclusion [23]. These cases are listed in the same descending order by the shear velocity * . Among the studied cases, the flow resistance in 3D-cu-s was the closest to the experi-  With the same flow rate, the roughness element arrangements had a great impact on the flow resistances. According to the values of k s /d and ∆U + , 2D-sq generated the largest flow resistance, followed by 3D-cu-c and 3D-cu-s, which corresponds to Volino et al.'s conclusion [23]. These cases are listed in the same descending order by the shear velocity u * . Among the studied cases, the flow resistance in 3D-cu-s was the closest to the experiment of Blanckaert et al., while the 2D-sq and 3D-cu-c profiles in Figure 7 were approximately coincident. It is not surprising that 3D-cu-s had the smallest flow resistance since its roughness elements were arranged in the sparsest way. In the region y + < 100, these velocity profiles deviated from the logarithmic law. Compared to the logarithmic law, the velocity-profile slopes for 2D-sq and 3D-cu-c in the y + < 100 region increased while that for 3D-cu-s decreased, indicating different vortex structures among these cases. For the case of 3D-cu-s, the horizontal vortices generated by the flow around the roughness elements were responsible for the nearly uniform velocity profile in the region y + < 100. Figure 8 displays the normalized Reynolds shear stress profiles and Figure 9 shows the normalized turbulent kinetic energy profiles, both of which were averaged horizontally over the computational domain. A sharp increase in the Reynolds shear stress in 2D-sq and 3D-cu-c was observed in the near-wall region, as the vertical recirculation flow in the cavities between the roughness elements produced a strong shear. A jump in the Reynolds stress profile was found at the top face of the roughness elements, above which the form drag due to the roughness element vanishment. Above the roughness element layer, the three-dimensional cases were roughly 20~30% lower than the two-dimensional case in the y + = 800~2600 range, which was noted by Volino et al., indicating the influence of the different roughness dimensions. There was a local peak for the Reynolds shear stress above the roughness elements, respectively, at y + = 237, y + = 86, and y + = 142 in 2D-sq, 3D-cu-c, and 3Dcu-s. The increasing and decreasing trends in the turbulent kinetic energy between the three cases were relatively similar, while their turbulent kinetic energy was basically dependent on the value of u * . However, the peak in 2D-sq was greater than in the two three-dimensional cases, indicating the difference in the flow structures among the three cases.

Characteristics of Time-Averaged Hydrodynamic Parameters
velocity-profile slopes for 2D-sq and 3D-cu-c in the 100 region increased while that for 3D-cu-s decreased, indicating different vortex structures among these cases. For the case of 3D-cu-s, the horizontal vortices generated by the flow around the roughness elements were responsible for the nearly uniform velocity profile in the region 100. Figure 8 displays the normalized Reynolds shear stress profiles and Figure 9 shows the normalized turbulent kinetic energy profiles, both of which were averaged horizontally over the computational domain. A sharp increase in the Reynolds shear stress in 2Dsq and 3D-cu-c was observed in the near-wall region, as the vertical recirculation flow in the cavities between the roughness elements produced a strong shear. A jump in the Reynolds stress profile was found at the top face of the roughness elements, above which the form drag due to the roughness element vanishment. Above the roughness element layer, the three-dimensional cases were roughly 20~30% lower than the two-dimensional case in the y + = 800~2600 range, which was noted by Volino et al., indicating the influence of the different roughness dimensions. There was a local peak for the Reynolds shear stress above the roughness elements, respectively, at y + = 237, y + = 86, and y + = 142 in 2D-sq, 3Dcu-c, and 3D-cu-s. The increasing and decreasing trends in the turbulent kinetic energy between the three cases were relatively similar, while their turbulent kinetic energy was basically dependent on the value of * . However, the peak in 2D-sq was greater than in the two three-dimensional cases, indicating the difference in the flow structures among the three cases.   ments were responsible for the nearly uniform velocity profile in the region 100. Figure 8 displays the normalized Reynolds shear stress profiles and Figure 9 shows the normalized turbulent kinetic energy profiles, both of which were averaged horizontally over the computational domain. A sharp increase in the Reynolds shear stress in 2Dsq and 3D-cu-c was observed in the near-wall region, as the vertical recirculation flow in the cavities between the roughness elements produced a strong shear. A jump in the Reynolds stress profile was found at the top face of the roughness elements, above which the form drag due to the roughness element vanishment. Above the roughness element layer, the three-dimensional cases were roughly 20~30% lower than the two-dimensional case in the y + = 800~2600 range, which was noted by Volino et al., indicating the influence of the different roughness dimensions. There was a local peak for the Reynolds shear stress above the roughness elements, respectively, at y + = 237, y + = 86, and y + = 142 in 2D-sq, 3Dcu-c, and 3D-cu-s. The increasing and decreasing trends in the turbulent kinetic energy between the three cases were relatively similar, while their turbulent kinetic energy was basically dependent on the value of * . However, the peak in 2D-sq was greater than in the two three-dimensional cases, indicating the difference in the flow structures among the three cases.   The influence of the rough wall to the time-averaged shear stress −u v could be partially reflected by its components u and v , while u rmse , v rmse , and w rmse are evaluations of the magnitude of the velocity fluctuations. Figure 10 displays the u rmse , v rmse , and w rmse profiles, normalized by u * . For the open-channel rough turbulence, u rmse + > w rmse + > v rmse + was observed [22,36]. The roughness elements exerted a greater impact on the velocity fluctuation in the near-wall region. Below y + = 100 in the roughness element layer, 2Dsq and 3D-cu-c had a similar u rmse + and v rmse + , while 2D-sq and 3D-cu-s had a similar w rmse + . This similarity is related to the overlapping of the normalized Reynolds stress in this region for 2D-sq and 3D-cu-c. In the outer region, the u rmse + differences between the cases were similar to that in the Reynolds stress profiles, with the three-dimensional cases around 12~22% lower while the v rmse + values were almost identical. This indicates that the differences in the Reynold stress might have mainly resulted from the differences in the longitudinal velocity fluctuation characteristics.
the magnitude of the velocity fluctuations. Figure 10 displays the urmse, vrmse, and wrmse p files, normalized by * . For the open-channel rough turbulence, urmse + > wrmse + > vrmse + w observed [22,36]. The roughness elements exerted a greater impact on the velocity fluc ation in the near-wall region. Below y + = 100 in the roughness element layer, 2D-sq a 3D-cu-c had a similar urmse + and vrmse + , while 2D-sq and 3D-cu-s had a similar wrmse + . T similarity is related to the overlapping of the normalized Reynolds stress in this region 2D-sq and 3D-cu-c. In the outer region, the urmse + differences between the cases were sim to that in the Reynolds stress profiles, with the three-dimensional cases around 12~2 lower while the vrmse + values were almost identical. This indicates that the differences the Reynold stress might have mainly resulted from the differences in the longitudi velocity fluctuation characteristics.  Volino et al. suggested that the complete obstruction over the channel width in tw dimensional cases forces flow detachment and acceleration over the bars and leads to larger ks/d and vertical range affected by the rough wall [23,37]. The complete blockage the lateral direction forces the flow to go over the roughness bars, offsetting the positio of the maximum time-averaged shear stress away from the wall [23]. This phenomeno was also observed in this study, as 2D-sq had the largest * , ks/d, and urmse + . Comparativel a horizontal flow around the roughness elements was allowed for the typical three-d mensional 3D-cu-s case, where uniform distributions of the horizontal velocity (both tim averaged values and their fluctuations) were found in the element roughness layer. Wit Volino et al. suggested that the complete obstruction over the channel width in twodimensional cases forces flow detachment and acceleration over the bars and leads to a larger k s /d and vertical range affected by the rough wall [23,37]. The complete blockage in the lateral direction forces the flow to go over the roughness bars, offsetting the position of the maximum time-averaged shear stress away from the wall [23]. This phenomenon was also observed in this study, as 2D-sq had the largest u * , k s /d, and u rmse + . Comparatively, a horizontal flow around the roughness elements was allowed for the typical three-dimensional 3D-cu-s case, where uniform distributions of the horizontal velocity (both time-averaged values and their fluctuations) were found in the element roughness layer. With regard to 3D-cu-c, the top surface drag of the roughness element forced its upstream flow to its two side cavities by going over the corner connections, although the chessboard-arranged roughness elements let no interstices between them. As a result, the flow in the roughness element layer was similar to 2D-sq, while the flow over the roughness element was similar to 3D-cu-s.
The overlapping of the u rmse + , v rmse + , and w rmse + profiles in the roughness element layer was another consequence of the roughness elements' blockage effect in each direction. In 2D-sq and 3D-cu-c, the arrangements of the roughness elements cut off the continuous longitudinal flow in the roughness element layer, which limited the longitudinal velocity fluctuations. In contrast, the arrangement of the roughness elements in 3D-cu-s allowed fewer limitations to both the longitudinal and lateral flow, resulting in a higher u rmse + than 2D-sq and 3D-cu-c and a similar w rmse + as 2D-sq in the roughness element layer, where flow in the lateral direction was also less limited.

Quadrant Analysis for Shear Stress
The scattered points on the u -v plane representing the type A and B locations each formed an elliptic cloud, and based on the number of the points that formed the cloud, the joint probability density function (PDF) of u and v can describe the probability distribution of fluctuations in each quadrant [38]. Here, 2D-sq is taken as an example and the joint PDF contours were calculated and are displayed in Figure 11. The major axes of the elliptic clouds and hyperbolas of the average shear stress are also plotted. On the u -v plane, the points of shear stress were mainly distributed in Q2 (36.38% in type A and 32.31% in type B) and Q4 (24.79% in type A and 36.05% in type B), with a stronger magnitude than the Q1 and Q3 fluctuations, indicating the dominance of the ejection and sweep event. The fluctuations in these two quadrants had positive contributions to the shear stress -u v , and for all three cases, the Q2 and Q4 fluctuations dominated and their total proportion in the type B points was higher than in the type A points, resulting in a higher average shear stress magnitude in the type B points. Table 3 provides the proportion of each quadrant in all three cases. In 2D-sq, the Q2 proportion was higher in the type A points, which is typical at the boundary layer edge from the front edge of the bars [38]. In the 3D cases, the Q4 proportions were higher in the type A points, as the arrangements of the staggeringly distributed cubes allowed velocity exchange between the type A points and the points beside them in the lateral direction with a higher-speed flow, which contributed to the Q4 fluctuations. For the type A and type B points, the two-dimensional case and three-dimensional cases had opposite top-dominating fluctuation events, which may explain the different profiles of the Reynolds shear stress between the two-dimensional case and three-dimensional cases. For the type A points in 2D-sq, the top-dominating event of the Q2 transported the larger shear stress upward and resulted in a smaller slope above the local peak point (Figure 8).
The fluctuations in Q2 and Q4 represent the exchange of energy in vertical directions and the accompanying mass transportation, such as the dislodgement of sediment and benthic macroinvertebrates. Previous research has revealed that Q2 fluctuations contribute to sediment dislodgement while Q4 fluctuations are related to mass transportations toward the bed. Q2 fluctuations contribute to particle dislodgement from the bed, which is partly related to the positive u' and thus the greater longitudinal velocity, and Q4 fluctuations maintain their suspension in the water column [8,9]. The spatial variation of the Q2 and Q4 events corresponds with previous research on microflow regimes around stream boulders, and is related to benthic macroinvertebrate distributions [39]. While shear stress is commonly regarded as a driving force of particle dislodgement from riverbeds, it may also benefit the transportation of particulate organic matter and the exchange of dissolved gases, which is another factor that relates to benthic macroinvertebrate distributions. Further analysis by decomposing shear stress components can build links to Stokes equations, converting shear stress components into lift force [9]. This builds a bridge between the flow shear stress and the actual force exerted on the particles, which allows an intensive understanding of sediment and benthic macroinvertebrate movements from a mechanical perspective.
each quadrant in all three cases. In 2D-sq, the Q2 proportion was higher in the type A points, which is typical at the boundary layer edge from the front edge of the bars [38]. In the 3D cases, the Q4 proportions were higher in the type A points, as the arrangements of the staggeringly distributed cubes allowed velocity exchange between the type A points and the points beside them in the lateral direction with a higher-speed flow, which contributed to the Q4 fluctuations. For the type A and type B points, the two-dimensional case and three-dimensional cases had opposite top-dominating fluctuation events, which may explain the different profiles of the Reynolds shear stress between the two-dimensional case and three-dimensional cases. For the type A points in 2D-sq, the top-dominating event of the Q2 transported the larger shear stress upward and resulted in a smaller slope above the local peak point (Figure 8).     Figure 12 displays the power spectrum of the shear stress temporal signal. The major part of the shear stress energy is contained in low-frequency fluctuations, and a −5/3-law can be observed in the power spectrum, which is similar to turbulent kinetic energy and scalar transportation [34,[40][41][42]. The frequency range of the −5/3-law is roughly from 8 Hz to 30 Hz. Taking the local time-averaged velocity as the characteristic velocity, the corresponding length scale range is about from 1 mm to 7 mm. The largest length scale has the same order of magnitude as the roughness element size and the smallest length scale has the same order of magnitude as the Kolmogorov length scale.

Spectrum Analysis for Shear Stress
allows an intensive understanding of sediment and benthic macroinvertebrate movements from a mechanical perspective. Figure 12 displays the power spectrum of the shear stress temporal signal. The major part of the shear stress energy is contained in low-frequency fluctuations, and a −5/3-law can be observed in the power spectrum, which is similar to turbulent kinetic energy and scalar transportation [34,[40][41][42]. The frequency range of the −5/3-law is roughly from 8 Hz to 30 Hz. Taking the local time-averaged velocity as the characteristic velocity, the corresponding length scale range is about from 1 mm to 7 mm. The largest length scale has the same order of magnitude as the roughness element size and the smallest length scale has the same order of magnitude as the Kolmogorov length scale. The power function F(f) added from low frequency to high frequency represents the accumulation of the shear stress energy from large scale vortices to small-scale vortices, and its maximum is the total energy of the shear stress temporal signal. In Figure 13, the frequency f is replaced with characteristic length scale , and the curves represent the accumulation of the shear stress signal energy from the largest scale to the smallest scale. The integral increases with the decrease in length scale, and it reaches the maximum at around 1 mm, corresponding to a frequency of around 30 Hz where the power spectral density is about two orders of magnitude smaller than its maximum. The power function F(f ) added from low frequency to high frequency represents the accumulation of the shear stress energy from large scale vortices to small-scale vortices, and its maximum is the total energy of the shear stress temporal signal. In Figure 13, the frequency f is replaced with characteristic length scale L v , and the curves represent the accumulation of the shear stress signal energy from the largest scale to the smallest scale. The integral increases with the decrease in length scale, and it reaches the maximum at around 1 mm, corresponding to a frequency of around 30 Hz where the power spectral density is about two orders of magnitude smaller than its maximum.  The total energy of the shear stress temporal signal was about 261.0~375.8% of the energy of time-averaged shear stress. According to the deduction of Equation (12), the theoretical value of this ratio should be 1 + u v − u v 2 /2 u v 2 , which is decided by the mean and the temporal fluctuation magnitude of the shear stress. The stronger the fluctuation (relative to the time-averaged shear stress), the larger the ratio. The power function revealed that the shear stress energy was mainly contained in fluctuations of 0~30 Hz, which corresponds to the energy cascade. For 3D-cu-s, the sparse distribution of the roughness elements resulted in low resistance to the flow, reflected by the lowest k s /d and u * . The rough wall in this three-dimensional case affected a relatively smaller range in the vertical direction ( Figure 10) [23]; therefore, it had a weaker influence on the large-scale fluctuations, leading to the smallest normalized shear stress energy in 3D-cu-s.

Spectrum Analysis for Shear Stress
In turbulent flows, large-size vortices contain the majority of the turbulent kinetic energy and are dominated by force from inertial range scales, exerting a pressure gradient force upon the particles in the fluid [43]. Small-size vortices dissipate the turbulent kinetic energy through viscosity, and the viscosity force is proportional to the torque upon the particles [44]. The available literature pointed out that for a finite-sized particle, as the Particle Reynolds number increases along with the particle size, the effect of small-size vortices decreases [43]. Experiments by Qureshi et al. revealed that the transport and rotation of particles are forced only by turbulent pressure fluctuations at scales larger than the particle [44].
By extending this argument to shear stress fluctuations, the accumulation of shear stress energy should be conducted from the largest length scale to the size of the object. The total effect of the flow fluctuations exerted upon the object decreases with the increase in its length scale. At the range of the benthic macroinvertebrate size (mm~cm), smaller benthic macroinvertebrates (~2 mm) are usually affected by all scales of vortices; therefore, small-scale flow fluctuations should not be neglected. For benthic macroinvertebrates and other particles of larger sizes (>2 mm), the effect of small vortices weakens with a scale increase; therefore, the effect from shear stress fluctuations is reduced. For larger objects in the flow, the effect of the shear stress fluctuation is even weaker, and finally the total effect of the shear stress is reduced to merely the effect of the time-averaged Reynolds stress which, together with the time-averaged flow velocity, has the main contribution.
When taking the difference of particle mass into consideration, fluctuations in turbulent flows may exert a stronger influence on small-size particle movement than expected. A previous experiment estimated that instantaneous velocity could be 40% higher than its time average, resulting in an increase of about 100% in lift force [9]. In this study, we evaluated the total energy of the shear stress signal, which could be 2~4 times the time-averaged part. This evaluation was based on the current flow field simulation, while the diversity of the turbulent flows and the complexity of fluid-solid coupled dynamics may result in different phenomena.

Conclusions
Both the mean value and the fluctuations of turbulent flow shear stress play critical roles in river ecological processes, but the fluctuation of shear stress has been seldom studied because of difficulties in direct measurements. In this study, we used numerical simulations to study shear stress characteristics over a natural rough riverbed, including the time-averaged and fluctuating characteristics. The roughness element method was used to mimic a natural rough bed. The arrangements of the roughness elements resulted in different degrees of blockage over the channel width by changing the flow pattern in the roughness element layer. Quadrant analysis revealed that ejections and sweep events were stronger behind the roughness element, which indicates a strong momentum exchange. This corresponds to the benthic macroinvertebrate distribution around boulders in the river. Spectrum analysis provides a method to quantitatively estimate the accumulated effect of instantaneous shear stress energy by vortex scales. With a decrease in the length scale, particles endure a larger influence from the fluctuating part of turbulence, and its magnitude could be about 2~4 times the time-averaged part of the flow. At scales of benthic macroinvertebrates, which is a key indicator of aquatic ecosystems, the effect of the instantaneous part of shear stress is not negligible and should be considered in ecological models.
The above fluctuation characteristics of shear stress are related to the dislodgement of sediment and benthic macroinvertebrates on the riverbed. The results and methods of this study will further support future intensive studies on river morphology studies and ecological processes, such as benthic macroinvertebrate drifting and bio-particle transport from a mechanical perspective.
Author Contributions: Data curation, J.W. and Z.L.; methodology, J.W. and Z.L.; supervision, Y.C.; writing-original draft, J.W.; writing-review and editing, Z.L. and M.L. All authors have read and agreed to the published version of the manuscript.