The E ﬀ ect of Habitat Structure Boulder Spacing on Near-Bed Shear Stress and Turbulent Events in a Gravel Bed Channel

: This study experimentally investigated the effect of boulder spacing and boulder submergence ratio on the near-bed shear stress in a single array of boulders in a gravel bed open channel flume. An acoustic Doppler velocimeter (ADV) was used to measure the instantaneous three-dimensional velocity components. Four methods of estimating near-bed shear stress were compared. The results suggested a signiﬁcant e ﬀ ect of boulder spacing and boulder submergence ratio on the near-bed shear stress estimations and their spatial distributions. It was found that at unsubmerged condition, the turbulent kinetic energy (TKE) and modiﬁed TKE methods can be used interchangeably to estimate the near-bed shear stress. At both submerged and unsubmerged conditions, the Reynolds method performed di ﬀ erently from the other point-methods. Moreover, a quadrant analysis was performed to examine the turbulent events and their contribution to the near-bed Reynolds shear stress with the e ﬀ ect of boulder spacing. Generally, the burst events (ejections and sweeps) were reduced in the presence of boulders. This study may improve the understanding of the e ﬀ ect of the boulder spacing and boulder submergence ratio on the near-bed shear stress estimations of stream restoration practices.


Introduction
Bed shear stress plays a determinant role in the incipient motion of sediment. The bed shear stress has been the focus of many studies in both laboratory flumes and rivers. The majority of the existing studies were conducted in the flumes due to the controlled flow and bed conditions. However, the field investigation on the measuring of the bed shear stress is relatively limited. The direct measurements on the bed shear stress were often conducted by using the shear plate [1][2][3][4][5]. The indirect measurement of the bed shear stress depends on the velocity in the inner layer, lower 20% of the flow depth [6]. The common approaches to calculate the bed shear stress include, but are not limited to, the reach-averaged bed shear [7], logarithmic law of the wall [8], drag force [9], Reynolds stress [10], and turbulent kinetic energy (TKE) methods [11]. In a laboratory study, the bed shear stress estimated from several methods were compared for a simple boundary layer and a complex flow [12]. It was found that the Reynolds and TKE methods showed the most appropriate results for a simple boundary layer and a complex flow, respectively [12]. In another study, Reynolds stress was calculated by using the double-averaged methodology based on Reynolds temporal and spatial averaging [13]. The study concluded that the maximum Reynolds stress could better represent the bed shear stress than the linearly extrapolated Reynolds stress profile [13]. In the field measurements, acoustic Doppler current

Scenario
Flow Rate (L/s)

Boulder Spacing
Reach-Averaged Flow Depth (m) Submergence Ratio S1- 60 Infinity 0.082 - Four series of experiments (hereafter called scenarios) were designed with boulder-to-boulder spacing (from center to center) as a varying parameter to represent different flow regimes. Covered boulder-to-boulder spacing was infinity (no boulder), 10D (large spacing), 5D (medium spacing), and 2D (small spacing). The large and medium boulder spacing referred to an isolated flow regime while the small boulder spacing corresponded to a wake-interference or even a skimming flow regime [38]. Table 1 lists four scenarios and the hydraulic parameters associated with each scenario. Each scenario was conducted at 60 L/s and 100 L/s flow rates. Boulders were unsubmerged at 60 L/s with a submergence ratio (H/D) between 0.73 and 0.78, while at 100 L/s, boulders were fully submerged with a submergence ratio between 1.25 and 1.32 (Table 1). Figure 2 illustrates the measurement points and position of boulders in each scenario in the measurement zone. For the boulder-scenarios, the boulders were placed along the centerline of the channel in a single straight array. First, the central boulder at point X4Y3 was placed; afterward, the other boulders were arranged relative to the central boulder. A total of 35 measurement points were selected covering an area from −4D to 4D (relative to the central boulder in the measurement zone) in the streamwise direction (x) and from −2D to 2D in the spanwise direction (y) of the central boulder. The distance of all the measurement points from the walls was greater than 10 cm to minimize the wall effects [28]. It is worth mentioning that 35 measurement points were measured for the no-boulder scenario; however, for the boulder-scenarios, some of the measurement points in the centerline of the boulders were blocked by the boulders. In addition, Figure 2 shows that the points immediately upstream and downstream of a boulder have an equal distance from the boulder; however, due to the irregularity in the boulder shape, these distances were not perfectly the same.

Data Filtering
Three-dimensional instantaneous velocities were collected using a Vectrino probe (Nortek ADV) at a sampling rate of 25 Hz for three minutes to ensure statistical significance. Aliased points were despiked using a phase-space threshold filter proposed by [49]. This filter can produce clean signals in highly turbulent flows [50]. The manufacturer (Nortek) recommends COR ≥ 70% and SNR ≥ 15 dB for reliable turbulence measurements. However, in the case of high SNR values-a valid condition in this study with an average SNR value of greater than 50 dB-COR values less than 70% can still provide reliable data [51]. As suggested by [52], COR values less than 70% in the turbulent flows are not necessarily an indicator of the low quality. They applied COR ≥ 50% and SNR ≥ 10 dB to The velocity measurements were taken in a relative depth (z/H) around 0.20-0.25, at the beginning of the outer region of flow [47]. Although measurements closer to the bed are preferable, the use of such a relative depth to capture the near-bed flow characteristics has been justified in a similar study [42]. Additionally, acoustic Doppler velocimeter (ADV) measurement about 3 cm from the bed to acquire higher levels of accuracy within the turbulent flows has been recommended [48], a threshold that led to a z/H around 0.20-0.30 in this study.

Data Filtering
Three-dimensional instantaneous velocities were collected using a Vectrino probe (Nortek ADV) at a sampling rate of 25 Hz for three minutes to ensure statistical significance. Aliased points were despiked using a phase-space threshold filter proposed by [49]. This filter can produce clean signals in highly turbulent flows [50]. The manufacturer (Nortek) recommends COR ≥ 70% and SNR ≥ 15 dB for reliable turbulence measurements. However, in the case of high SNR values-a valid condition in this study with an average SNR value of greater than 50 dB-COR values less than 70% can still provide reliable data [51]. As suggested by [52], COR values less than 70% in the turbulent flows are not necessarily an indicator of the low quality. They applied COR ≥ 50% and SNR ≥ 10 dB to investigate the turbulent characteristics around a cluster microform [52]. It also has been shown that for the flow over a rough bed, if at least 70% of data remain after applying a filtering scheme, the Reynolds stress could be measured with a minimum COR of 40% [53]. Following these suggestions, a filtering scheme with COR ≤ 65% and SNR ≤ 20 dB was applied to filter out poor quality data. For the points in location X5Y3 (marked in Figure 2a) at 100 L/s flow rate, a COR of 55% was applied due to the higher turbulence in this point and a significant loss of data in case of applying the main filtering scheme. Such criteria resulted in the remaining 77.1% of data indicating an acceptable portion of data for further analysis. Another possible source of error in ADV data is the Doppler noise, which may affect the turbulence parameters [54]. When the noise level is low, it is expected that the velocity spectra follow the Kolmogorov's −5/3 law by exhibiting a −5/3 slope in the inertial subrange of the velocity spectra [54]. Therefore, as an additional check (also suggested by [50] for highly turbulent flows), the velocity spectra of the measured time series were visually inspected to remove noisy time series with a flat or negative slope in the inertial subrange. None of the points were rejected based on this criterion.

Calculation of Bed Shear Stress
The first main method to calculate the bed shear stress is the reach-averaged method as below [12]: where τ 0 is bed shear stress, ρ is water density, g is gravitational acceleration, and R is channel hydraulic radius. Although this method gives a rough estimation of bed shear stress, its performance is questionable, specifically in the presence of boulders, because it is not able to estimate the local variation of the bed shear stress [12]. Bed shear stress can be estimated directly through the Reynolds shear stress [12]: where u and w are the instantaneous velocity fluctuations of the streamwise and vertical components, respectively. The bar sign over u w indicates a time-averaged value. This method is sensitive to the deviations from a two-dimensional uniform flow [12]. Another method is the turbulence kinetic energy (TKE) method as below [12]: where v is the instantaneous velocity fluctuation of the spanwise component and C 1 is a constant equal to 0.19 [11]. The modified TKE method only uses the vertical component of the velocity fluctuations due to the smaller noise of the measurement instrument in this direction [11]. Therefore, it can be written as below [12]: where C 2 is a constant equal to 0.9 [11]. Equations (1)-(4) were used to calculate the near-bed shear stress for all the conducted scenarios in this study. Equation (1) only requires the reach hydraulic data while Equations (2)-(4) use the measured data from a single point and are referred to as point-methods in this study.

Quadrant Analysis
In a quadrant analysis, the velocity fluctuations are decomposed to four quadrants based on their signs [38]: To remove insignificant events, a hole region in the (u , w ) plane can be defined to exclude points within this region. Hole size, H, is a parameter to define the expansion of the hole region as below [55]: Clearly, a larger hole size means the remaining of more significant events. Many studies have selected a hole size greater than one to distinguish significant events [36,[56][57][58][59]. Therefore, in this study, a hole size of 2 was selected to distinguish events with a higher magnitude. Points X3Y3 and X5Y3, the common points in all the scenarios (marked in Figure 2a), were selected for a quadrant analysis. The joint frequency distribution of the normalized velocity fluctuations was found for each point. To find the contribution of each quadrant to the Reynolds shear stress, the stress fraction parameter of each quadrant, S f i.H was found as below [38]: where i shows the quadrant number (i = 1, 2, 3, 4), T is the total duration of the measurement and I is a conditional function to detect events in each quadrant for a given hole number as below [38]: The total portion of each turbulent event, P i,H , as a representative of the frequency of turbulent events can be determined as below [38]:

Near-Bed Shear Stress
Equations (1)-(4) were used to calculate the reach-averaged near-bed shear stress for each scenario. Equation (1) gives a single value as the reach-averaged shear stress, but for the point-methods, the average of shear stresses at all the measurement points was calculated as the reach-averaged near-bed shear stress. Figure 3 shows the reach-averaged near-bed shear stress for each method and scenario. The near-bed shear stress estimations using the reach-averaged method (Equation (1)) before adding boulders were 3.46 and 5.63 N/m 2 for the S1-60 and the S1-100 scenarios, respectively. These estimations were higher than the shear stress estimations using the point-methods for no-boulder scenarios. The maximum estimation of the average near-bed shear stress was significantly lower than the calculated critical bed shear stress (12.38 N/m 2 ). This was expected because no sediment movement was observed in all scenarios. For the no-boulder scenarios (S1-60 and S1-100), the TKE and modified TKE methods resulted in a close estimation for the average near-bed shear stress (a maximum difference of 10% for the S1-60) while estimations using the Reynolds method were slightly higher than the TKE and modified TKE methods. Close estimations using these point-methods were reported by [12]. For scenarios with the largest boulder spacing at 60 L/s (S2-60), the TKE method resulted in a higher bed shear stress while the Reynolds and modified TKE resulted in close estimations. At 100 L/s scenario (S2-100), the Reynolds, TKE, and modified TKE methods showed very similar performance (with a maximum 20% difference between the modified TKE and TKE methods). At 60 L/s, for both medium (S3-60) and small boulder spacing (S4-60) scenarios, both TKE and modified TKE methods showed a very similar performance while the Reynolds method showed a significantly lower near-bed shear stress. The average near-bed shear stress using the Reynolds method was about 200% and 70% less than estimations of other point-methods for the S3-60 and S4-60, respectively. At 100 L/s, for both S3-100 and S4-100, the modified TKE method led to the highest estimates, while the Reynolds and TKE methods showed similar results. Generally, it can be said that at submerged condition (100 L/s), the difference between estimations using the point-methods was not significant for each scenario, although the modified TKE led to slightly higher estimations for boulder-scenarios. At unsubmerged condition (60 L/s), for both no-boulder and large spacing scenarios, the difference between the point-methods was not significant; however, for the medium and small spacing scenarios, the near-bed shear estimations using the Reynolds method was significantly smaller.   Figure 5 shows plots for comparing estimations using the Reynolds shear stress method versus the TKE and modified TKE methods. At 60 L/s and 100 L/s, before adding boulders (S1-60 and S1-100), both TKE and modified TKE methods showed a reasonable agreement with the Reynolds method (At 60 L/s, RMSE = 0.63 and 0.73 N/m 2 for the TKE and modified TKE, respectively. At 100 L/s, RMSE = 0.48 and 0.43 N/m 2 for the TKE and modified TKE, respectively). It was reported that the modified TKE and Reynolds method have a relatively good agreement over a Plexiglas and sand bed [12]. However, estimations using the TKE and modified TKE methods were mostly lower than the Reynolds method, but at a lower range of shear stresses, the values were close to 1:1 line. This behavior in a range of high shear stress is in agreement with [12]; however, for the lower shear stresses, similar studies found systematically higher values from the TKE and modified TKE methods [12,40]. For the large boulder spacing at 60 L/s (S2-60), the agreement between the Reynolds and both TKE and modified TKE methods slightly decreased (RMSE = 0.87 and 0.95 N/m 2 for the TKE and The effect of adding boulders and boulder spacing can also be interpreted in Figure 3. The reach-averaged method is directly related to the average flow depth (Equation (1)), and a higher average flow depth resulted in a higher bed shear stress. After adding boulders, the average flow depth slightly increased, and the reach-averaged bed shear stresses slightly increased too. However, for the S4-60 and S4-100 scenarios, the flow depth, and subsequently, the reach-averaged bed shear stress slightly decreased. For the Reynolds method at 60 L/s scenarios, it can be seen that the bed shear stress gradually decreased by decreasing boulder spacing-except for the S4-60 scenario, in which the bed shear stress slightly increased relative to the previous scenario at S3-60. The near-bed shear stress variation for different boulder spacing at 60 L/s using the TKE and modified TKE methods was not similar to the Reynolds method. For the TKE method, the near-bed shear stress slightly increased for the large boulder spacing scenario (S2-60) and remained almost constant for the medium spacing scenario (S3-60); then, it reached the minimum for the small spacing scenario (S4-60). For the modified TKE method, the highest near-bed shear stress belonged to the medium spacing scenario (S3-60), while the other scenarios had a small difference with each other (less than 15%). At 100 L/s, the small spacing scenario (S4-100) experienced the highest near-bed shear stress resulting from all estimation methods. For the Reynolds and TKE methods, the near-bed shear stress decreased for the large spacing scenario (S2-100) and slightly increased for the medium spacing scenario (S3-100); then, it sharply increased for the scenario with the smallest boulder spacing (S4-100). For the modified TKE method, the near-bed shear stress gradually increased after adding boulders and then decreasing boulder spacing; the sharpest increase occurred for the small spacing scenario (S4-100). This trend is in agreement with the findings of [38] in which the total shear stress increased gradually by increasing the number of boulders. In the following sections, a more detailed discussion on the variation of near-bed shear stress for different boulder spacing and estimation methods can be found.

Relative Performance of the Reynolds, TKE, and Modified TKE Methods
To compare the estimates of Reynolds, TKE, and modified TKE methods in a more compelling way, the near-bed shear stresses from each method for all the points were plotted against each other. Figure 4 shows the near-bed shear stress estimated using the TKE and modified TKE methods for all scenarios. Here, the goal was to understand the similarity between the near-bed shear stress estimated from the point-methods. Therefore, the root mean square error (RMSE) was selected to shed light on the similarity of the results calculated by different methods. As shown in Figure 4a, at 60 L/s scenarios, the RMSE for all the scenarios varied between 0.36 and 0.84 N/m 2 indicating a reasonable similarity between results even after adding boulders and then decreasing the boulder spacing. At 100 L/s (Figure 4b), before adding boulders (S1-100), a low RMSE = 0.18 N/m 2 indicated a good agreement between the TKE and modified TKE methods. After adding boulders, the RMSE significantly increased, and this increase was intensified after decreasing the boulder spacing. Generally, the modified TKE method resulted in higher values of the near-bed shear stress in comparison with the TKE method for the boulder-scenarios (S2-100, S3-100, and S4-100). This difference can be related to the values of C 1 and C 2 constants that still are required to be found for the natural streams [11,12]. Figure 5 shows plots for comparing estimations using the Reynolds shear stress method versus the TKE and modified TKE methods. At 60 L/s and 100 L/s, before adding boulders (S1-60 and S1-100), both TKE and modified TKE methods showed a reasonable agreement with the Reynolds method (At 60 L/s, RMSE = 0.63 and 0.73 N/m 2 for the TKE and modified TKE, respectively. At 100 L/s, RMSE = 0.48 and 0.43 N/m 2 for the TKE and modified TKE, respectively). It was reported that the modified TKE and Reynolds method have a relatively good agreement over a Plexiglas and sand bed [12]. However, estimations using the TKE and modified TKE methods were mostly lower than the Reynolds method, but at a lower range of shear stresses, the values were close to 1:1 line. This behavior in a range of high shear stress is in agreement with [12]; however, for the lower shear stresses, similar studies found systematically higher values from the TKE and modified TKE methods [12,40]. For the large boulder spacing at 60 L/s (S2-60), the agreement between the Reynolds and both TKE and modified TKE methods slightly decreased (RMSE = 0.87 and 0.95 N/m 2 for the TKE and modified TKE, respectively). By decreasing boulder spacing (S3-60 and S4-60), the RMSE increased to higher values indicating a lower similarity between the methods. A reason for the poor performance can be the presence of negative values in the Reynolds method, while values using two other methods are always positive. After adding boulders, it can be found that both TKE and modified TKE methods led to higher estimations than the Reynolds method (except for the modified TKE method at the large boulder spacing). This difference became more significant for the medium and small boulder spacing scenarios.
At 100 L/s for the boulder-scenarios, the RMSE of the Reynolds and TKE methods varied between 1.16 and 1.98 N/m 2 , and the RMSE of the Reynolds and modified TKE methods varied between 1.30 to 2.67 N/m 2 . The RMSE values for both TKE and modified TKE methods reached their maximum (lower similarity) for the small boulder spacing (S4-100). However, some extreme values resulted from the TKE and modified TKE methods might be the reason behind higher RMSE values. After adding boulders, the TKE method estimates generally were smaller than the Reynolds method, while the modified TKE method usually resulted in higher near-bed shear stress than the Reynolds method.

Spatial Distribution of Near-Bed Shear Stress
The spatial distribution of the near-bed shear stress is also of great importance, especially in the presence of boulders that significantly alter local flow fields. On this account, contour maps of estimated local near-bed shear stress using different methods were plotted. Figure 6 shows contour plots of the estimated near-bed shear stress using the Reynolds, TKE, and modified TKE methods at unsubmerged condition (60 L/s). Before adding boulders (S1-60) for all methods, a uniform spatial distribution of the near-bed shear stress was expected. However, as shown in Figure 6, slightly higher near-bed shear stresses can be seen on the left side of the measurement zone before adding boulders.

Spatial Distribution of Near-Bed Shear Stress
The spatial distribution of the near-bed shear stress is also of great importance, especially in the presence of boulders that significantly alter local flow fields. On this account, contour maps of estimated local near-bed shear stress using different methods were plotted. Figure 6 shows contour plots of the estimated near-bed shear stress using the Reynolds, TKE, and modified TKE methods at unsubmerged condition (60 L/s). Before adding boulders (S1-60) for all methods, a uniform spatial distribution of the near-bed shear stress was expected. However, as shown in Figure 6, slightly higher near-bed shear stresses can be seen on the left side of the measurement zone before adding boulders. This non-uniform distribution can be attributed to the presence of microforms, placed slightly higher than the mean bed level, on some points of the gravel bed in these series of experiments. For the Reynolds method, after adding boulders, it can be observed that the near-bed shear stress generally decreased in all measurement points. For the large spacing scenario (S2-60), a zone with a highly decreased near-bed shear stress can be observed extended to a distance of 2D downstream of the boulder. For the medium and large boulder spacing scenarios (S3-60 and S4-60), highly decreased shear stresses can be observed along the boulder centerline and extended to the sides of the boulder. It is known that the presence of an unsubmerged boulder within the flow generates intense reverse and irregular flow in the extended wake region [30]. Here, negative shear stresses downstream of the boulder can be attributed to this reverse flow [60]. Additionally, at the smaller boulder spacing with a wake-interference or skimming flow regime, the flow deceleration becomes more pronounced between the boulder [38]. This can be the reason for the low shear stress zone along the boulder centerline for scenarios with a smaller boulder spacing (S3-60 and S4-60) [38]. The spatial distribution of the near-bed shear stress using the TKE and modified TKE methods were similar to each other. The reason for this similarity can be a less significant streamwise velocity due to the reverse flow that led to a more important role of w in the calculation of shear stress using the TKE and modified TKE shear stress. Therefore, in this condition, the behavior of the TKE method became more similar to the modified TKE method that only depends on w . However, the TKE and modified TKE spatial distribution deviated from the Reynolds method. In contrast to the Reynolds method, there was a slight increase in the near bed-bed shear stress downstream of the boulders (except downstream of one of the boulders in the S4-60). This increase became more significant for the medium boulder spacing. On the sides of the boulder, the near-bed shear stress did not vary significantly in comparison with the no-boulder scenario.
At 100 L/s (Figure 7), the spatial distribution of the near-bed shear stress using the point-methods was similar to each other, but it was different from 60 L/s. In the large boulder spacing scenario (S2-100), the near-bed shear stress decreased upstream of the boulder, but a significant increase can be seen downstream of the boulder (extended to around 2D) for all the methods. Studies have reported an increased Reynolds shear stress and TKE downstream of a submerged boulder near the bed due to the large-scale vortices in the wake zone and the presence of vertical eddies due to the downwash flow [21,31,61,62]. This increased near-bed shear stress decreased with increasing downstream x/D distance. Similarly, it has been reported this high shear stress zone extended up to the reattachment point, and then the shear stress started to decrease as undisturbed turbulent intensities recovered [31,62]. For the medium boulder spacing scenarios (S3-100), again an increased shear stress zone can be observed downstream of the central boulder for all the methods; however, for the Reynolds method a decrease in the near-bed shear stress observed at X1Y3 at x/D = −4D, downstream of one of the boulders (this boulder has not been shown in the contour plot). Additionally, the length of the high shear stress zone decreased to approximately 1D. Upstream of the boulder, the near-bed shear stress using the Reynolds method decreased, while increased shear stress can be seen upstream of the boulders for the TKE and modified TKE methods. For the small boulder spacing (S4-100), high shear stress zones became more pronounced in the sides of the boulder in this scenario. Increased near-bed shear stresses can also be seen between boulders probably as a result of the stable vortices between boulders in a wake-interference or skimming flow regime [63]. However, for the Reynolds method, a decreased shear stress zone was observed upstream of the central boulder (X3Y3 at x/D = −1D). The reason for this alternate variation of the Reynolds method estimations between boulders for the medium and small boulder spacing remained unclear and required additional experiments. For all scenarios at 100 L/s, it can be seen that the modified TKE method resulted in a higher magnitude of near-bed shear stress, as already shown in Figure 3b. The main reason for this can be attributed to the uncertain values of C 1 and C 2 , as mentioned earlier.

Near-Bed Turbulent Events
Quadrant analysis was conducted for the common points in all scenarios, X3Y3 and X5Y3. Figures 8 and 9 show the joint frequency distribution of normalized ' and ' at X3Y3 and X5Y3, respectively. ' and ' were normalized by the reach-averaged velocity of each scenario. Tables 2 and 3 complement the results of the quadrant analysis by listing the stress fraction parameter ( ) and the frequency parameter ( ) for each quadrant at X3Y3 and X5Y3, respectively. Before adding boulders at both points and flow rates, as expected ejection and sweep events were dominant near stress and frequency. This is in agreement with the findings of [38] in which they found ejections and sweeps dominant for 10D and 5D spacing, while for 2D spacing each quadrant showed an almost equal frequency upstream of a fully submerged boulder [38]. This can be attributed to the stable vortices and a less pronounced downwash flow at the wake-interference or skimming flow associated with the small boulder spacing.

Near-Bed Turbulent Events
Quadrant analysis was conducted for the common points in all scenarios, X3Y3 and X5Y3. Figures 8 and 9 show the joint frequency distribution of normalized u and w at X3Y3 and X5Y3, respectively. u and w were normalized by the reach-averaged velocity of each scenario. Tables 2 and 3 complement the results of the quadrant analysis by listing the stress fraction parameter (S f i ) and the frequency parameter (P i ) for each quadrant at X3Y3 and X5Y3, respectively. Before adding boulders at both points and flow rates, as expected ejection and sweep events were dominant near the bed, and this was in agreement with the previous studies [41,64]. At both flow rates at X3Y3, after adding and then decreasing the boulder spacing, the shear stress fraction of ejections and sweeps increased while the frequency of them gradually decreased. However, sweeps and ejections remained dominant for the large and medium boulder spacing. As can be seen in Figure 8 and Table 2, for the small boulder spacing, all four quadrants showed almost an equal contribution to the Reynolds shear stress and frequency. This is in agreement with the findings of [38] in which they found ejections and sweeps dominant for 10D and 5D spacing, while for 2D spacing each quadrant showed an almost equal frequency upstream of a fully submerged boulder [38]. This can be attributed to the stable vortices and a less pronounced downwash flow at the wake-interference or skimming flow associated with the small boulder spacing.
very small. However, as shown in Table 3, between ejection and sweep events, sweeps showed significantly higher and . This was not in agreement with the findings of [38] specifically for 5D and 2D boulder spacing [38]. Besides, and showed a consistent trend that means kept the same order of dominance of (e.g., when for ejections and sweeps was the highest, was also the highest for ejections and sweeps). The turbulent events have a considerable effect on sediment transport and resuspension. Predominant burst events (ejections and sweeps) may affect the sediment resuspension and transport around the boulder [39][40][41]. Specifically, downstream of the boulder, sweeps may remain effective for the bedload transport at the medium and large boulder spacing, while less significant ejections may decrease the sediment entrainment and the fine sediment suspension in these zones. Additionally, adding boulders and reduction of ejection and sweep events may provide resting zones and a path to avoid high ejection zones for the fish [43]. Such areas can be found downstream of the boulders at the medium and large boulder spacing where the burst events reduced in comparison with the no-boulder condition. However, more investigation is necessary to find the role of turbulent events in sediment transport and habitat preference for different boulder spacing.      At 60 L/s at X5Y3 (Figure 9a), for the large boulder spacing (S2-60), inward and outward interactions were dominant. At the medium boulder spacing (S3-60), outward interactions and sweep events showed higher P i and S f i . Again, at the small boulder spacing (S4-60), outward and inward interactions became dominant. Although 60 L/s indicates an unsubmerged boulder condition, the pattern of changes in the turbulent events was completely in agreement with the results of [38] for the different spacing of a fully submerged boulder. At 100 L/s at X5Y3 (Figure 9b), ejection and sweep events remained dominant in all scenarios while the share of inward and outward interactions was very small. However, as shown in Table 3, between ejection and sweep events, sweeps showed significantly higher P i and S f i . This was not in agreement with the findings of [38] specifically for 5D and 2D boulder spacing [38]. Besides, P i and S f i showed a consistent trend that means P i kept the same order of dominance of S f i (e.g., when S f i for ejections and sweeps was the highest, P i was also the highest for ejections and sweeps). The turbulent events have a considerable effect on sediment transport and resuspension. Predominant burst events (ejections and sweeps) may affect the sediment resuspension and transport around the boulder [39][40][41]. Specifically, downstream of the boulder, sweeps may remain effective for the bedload transport at the medium and large boulder spacing, while less significant ejections may decrease the sediment entrainment and the fine sediment suspension in these zones. Additionally, adding boulders and reduction of ejection and sweep events may provide resting zones and a path to avoid high ejection zones for the fish [43]. Such areas can be found downstream of the boulders at the medium and large boulder spacing where the burst events reduced in comparison with the no-boulder condition. However, more investigation is necessary to find the role of turbulent events in sediment transport and habitat preference for different boulder spacing.

Conclusions
This study investigated the effects of boulder spacing and boulder submergence ratio on the near-bed shear stress estimations and their distributions using laboratory experiments. The following conclusions can be drawn from this study:

1.
For all scenarios, the reach-averaged method led to a significantly higher reach-averaged bed shear stress. For the unsubmerged condition, the Reynolds method resulted in a significantly lower near-bed shear stress between the point-methods, while, at submerged condition, all the point-methods showed very similar results.

2.
At unsubmerged condition, the effect of the boulder spacing on the variation of near-bed shear stress estimated from the Reynolds method was different from the TKE and modified TKE methods, while, at submerged condition, all of the point-methods showed a similar trend.

3.
At submerged condition, the densest boulder spacing led to the highest near-bed shear stress for all point-methods. However, for the unsubmerged condition, maximum near-bed shear stress varied for different methods and boulder spacing.

4.
At unsubmerged condition, the TKE and modified TKE methods can be used interchangeably for estimation of the near-bed shear stress in the presence of boulders; however, applying appropriate C 1 and C 2 coefficients is required to obtain more reliable results. For a comprehensive comparison of the Reynolds method with two other point-methods, more measurements are necessary, especially at unsubmerged condition. 5.
At unsubmerged condition, denser boulder spacing led to a more uniform contribution of turbulent events to the Reynolds shear stress. At submerged condition, decreased ejection events downstream of boulders in the large and medium boulder spacing may reduce the sediment entrainment and suspension.
This work is based on two boulder submergence ratios and single-point measurements near the bed and placing boulders along the centerline of the channel in a single straight array. To derive general conclusions, further measurements closer to the bed with a denser spatial grid between boulders are necessary. Additionally, a study of the bedload and the suspended load sediment transport within different boulder setups is recommended as the next step. Funding: This research received no external funding. However, financial support for the first author was provided by Clarkson University.