A Comparative Analysis of Logarithmic Law, Reynolds Stress, and Turbulence Kinetic Energy Methods for Bed Shear Stress Estimation in Complex Open‐Channel Flow

Numerous indirect methods for estimating the bed shear stress using velocity or turbulent stress profiles have been suggested in previous research. Although these methods have proven effective for simple boundary‐layer‐type flows, their efficacy in complex flow scenarios remains largely unexplored. This study aimed to evaluate the predictive capabilities of three popular indirect bed‐shear‐stress estimation methods—the logarithmic law, Reynolds shear stress (RSS), and turbulence kinetic energy (TKE) techniques—in a complex flow environment involving an obstacle in an open channel, producing massive flow separation and unsteady vortex shedding. To circumvent the difficulties of direct bed shear‐stress measurements, the reference bed shear stress was obtained from a high‐resolution wall‐resolving large‐eddy simulation (LES) data set. The key findings of this study are as follows: First, the logarithmic law and TKE methods were effective only in regions where the streamlines were almost parallel to the primary flow direction. Second, the ratio of the bed shear stress to TKE varied significantly in space in complex‐flow regions, rendering TKE methods virtually ineffective in these areas. Third, the RSS methods successfully reproduced the LES‐computed bed shear stress distributions, both qualitatively and quantitatively. Fourth, the accuracy of the RSS methods was influenced by two critical factors: (a) the incorporation of the transverse RSS component in the RSS extrapolation and (b) the selection of extrapolation techniques. Finally, this study recommends the use of RSS methods employing two‐point extrapolation for bed shear‐stress estimation in complex flows.


Introduction
The interaction between flow and sediment transport in rivers causes morphological changes in the bed, riverbank migration, and local scouring around hydraulic structures.Among the various factors, bed shear stress is widely recognized as the most critical flow parameter for understanding flow-sediment interactions.Consequently, in the fields of fluvial geomorphology and river engineering, considerable research has been dedicated to developing effective methodologies for the direct measurement or approximation of bed shear stress in openchannel flows.
A straightforward way to estimate the bed shear stress is to calculate the reach-average bed shear stress (Babaeyan-Koopaei et al., 2002;Biron et al., 2004;Henderson, 1966;Smart, 1999;Yang et al., 2015) as follows: where τ 0 is the bed shear stress, ρ is the density of water, g is the gravitational acceleration, R is the hydraulic radius, and S f is the energy slope due to friction.This approach is based on the equilibrium between gravitational force and boundary shear stress in a uniform open-channel flow.Therefore, it is unsuitable for the estimation of local bed stress (Biron et al., 2004;Smart, 1999), given the potential differences between the local and average stresses under nonuniform flow conditions.Several prior studies have attempted to formulate analytical expressions aimed at estimating bed shear stress under nonuniform flow conditions in rectangular or prismatic compound channels (Dey & Lambert, 2005;Guo & Julien, 2005;Liu et al., 2013;Tang & Knight, 2008;Yang & Lim, 1997).However, a universally applicable analytical method for all channel types is yet to be developed.
Numerous efforts have been made in the fields of oceanography (Conley & Griffin, 2004;Gust, 1988;Pujara & Liu, 2014;Rankin & Hires, 2000) and hydraulics (Cardoso et al., 1991;Park et al., 2016;Pasternack et al., 2007) to directly measure the local bed shear stress in open-channel settings, particularly by employing instruments such as hot films, shear plates, and torque sensors.Although the direct measurement approach enables precise measurements of the wall shear stress, it requires careful calibration and installation of instruments.Moreover, the measurement results are highly sensitive to the position and orientation of the sensor, its proximity to roughness elements, and temperature variations (Graham et al., 1992;Gust, 1988;Pope et al., 2006;Rankin & Hires, 2000).Consequently, direct measurement techniques have primarily been used to determine the shear stress at a few points in carefully controlled laboratory settings.However, performing such measurements in field environments may be impractical.The challenge of directly measuring bed shear stress with reliable accurately has led to the development of indirect methods.These methods aim to infer bed shear stress by correlating time-averaged velocity or turbulent stress profiles with bed shear stress values.The basis for these methods is that bed shear stress can be estimated using near-bed turbulence statistics, as they demonstrate a strong correlation.The advantage of these methods is that they enable the identification of bed shear stress distributions across large areas with significantly less effort compared to direct methods.
In the literature, three different indirect methods are identified: the logarithmic law method, the Reynolds shear stress (RSS) method, and the turbulence kinetic energy (TKE) method.Among these, the logarithmic law method is the most well-known, wherein the local bed shear is determined using the logarithmic law of the wall (Cardoso et al., 1989;Jackson, 1981;Nezu & Rodi, 1986;Yang et al., 2015): where U is the mean velocity, u* is shear velocity (= ̅̅̅̅̅̅̅̅ ̅ τ 0 /ρ √ ), κ is von Karman's constant (∼0.4), z is the height above the bed, and ν is kinematic viscosity.B is a constant that depends on relative roughness (Cardoso et al., 1989;Schlichting & Gersten, 2000;Wilcock, 1996;Yang et al., 2015).According to Pope (2000) and Yang et al. (2004), the value of constant B in a plain channel flow is ∼5.Theoretical and experimental studies have confirmed the effectiveness of this method for flow in a straight channel without obstruction.However, this method is ineffective for complex flows (Biron et al., 2004).
Alternatively, the local bed shear stress can be estimated by extrapolating the RSS to the bed (Graf & Song, 1995;Song & Chiew, 2001;Song et al., 1994): where u′ and w′ are the velocity fluctuations in the streamwise and vertical directions, respectively, and 〈〉denotes time averaging.This methodology, hereafter referred to as the RSS method, is based on the hypothesis that the RSS dominates the viscous stress near the bed and that the peak RSS value, which typically occurs near the bed, serves as a reliable estimate of the local bed shear stress.This has been theoretically demonstrated for plain channel flows and zero-pressure-gradient boundary-layer flows.Biron et al. (2004), Dey and Barbhuiya (2005), Duan (2009), Duan et al. (2009), Kuhnle et al. (2008), andJeon et al. (2018) used this method to estimate the bed shear stress in complex open-channel flows involving obstacles.
Another category of indirect methods that is particularly popular in oceanography is the TKE method.This technique postulates a linear correlation between the bed shear stress and TKE near the bed (Huthnance et al., 2002;Kim et al., 2000;Soulsby, 1983;Stapleton & Huntley, 1995;Williams et al., 1999) and calculates the shear stress using the following relation: where k is the TKE and C 1 is an empirical constant of approximately 0.19 (Soulsby, 1983;Stapleton & Huntley, 1995).This method has been used in oceanographic and fluvial studies (Biron et al., 2004;Hopkinson & Wynn-Thompson, 2012;MacVicar & Roy, 2007, 2011;Noss et al., 2010;Pope et al., 2006;Stapleton & Huntley, 1995;Thompson et al., 2003).The main limitation of the TKE method is that Equation 4 was not originally established based on rigorous theoretical analysis, but rather derived from empirical observations.For example, a constant bed shear stress-to-TKE ratio near the bed has been observed in field measurements of tidal boundary layer flows (Kim et al., 2000;Soulsby, 1983).

Water Resources Research
10.1029/2023WR036153 Although indirect methods have been applied to a variety of open-channel flow applications in the past, their primary application has been straight open-channel flows without obstructions.In these cases, the logarithmic law can provide accurate bed shear-stress values.However, the predictive capability of indirect methods for complex flows that often involve obstacles, particularly in comparison to direct bed shear stress measurements, remains largely unexplored.This is presumably because accurate shear stress measurements are difficult to obtain under such conditions.Nevertheless, several studies have evaluated the performance of different indirect methods, frequently utilizing velocity data obtained through Acoustic Doppler Velocimetry (ADV) (Biron et al., 2004;MacVicar & Roy, 2007;Pope et al., 2006;Thompson et al., 2003).Biron et al. (2004) compared the logarithmic law, RSS, and TKE methods for a straight open-channel flow with and without two past deflectors on the sidewalls in a laboratory flume under both fixed and mobile bed conditions.They showed that all three methods produced similar results in the case of straight flow over a fixed bed, particularly when velocity sampling was performed at approximately 10% of the depth above the bed.In a mobile-bed scenario with deflectors, the authors observed that the bed shear stress estimated using the TKE method exhibited a stronger correlation with the measured bed topography than the RSS method.Thompson et al. (2003) compared various indirect methods using velocity measurement data.MacVicar and Roy (2007) calculated the local bed shear stress in a gravel-bed setting using different indirect methods, and compared their results with those obtained using the energy slope equation (Equation 1).Despite the research, much remains to be done to clarify the capabilities and limitations of the indirect methods.It is important to determine whether indirect methods can effectively replicate bed shear-stress patterns under complex flow settings involving complex in-stream obstacles or irregular bathymetry.
The objective of this research was to evaluate the capability of various indirect methods for estimating bed shear stress in complex flow environments.To this end, this study considered the flow past a rectangular obstacle attached to the side of a laboratory flume, which was studied both experimentally and numerically by Jeon et al. (2018) and Kang et al. (2023Kang et al. ( , 2021)), respectively.Jeon et al. (2018) performed velocity measurements around and in the wake of an obstacle using an ADV.Later, Kang et al. (2021) performed a high-resolution largeeddy simulation (LES) using approximately 160 million grid points with near-bed spacings of less than one wall unit.These studies revealed that strongly three-dimensional flow patterns and complex flow features develop around and downstream of obstacles.These include horseshoe vortex systems, recirculation zones, energetic vortex shedding, turbulent shear layers, and secondary flows.The LES-computed flow statistics (Kang et al., 2021(Kang et al., , 2023) ) were validated against ADV measurements (Jeon et al., 2018) and flow visualization results (Kang et al., 2023), demonstrating satisfactory agreement.This represents a unique data set that offers both accurate flow field and bed shear stress data simultaneously, which is essential for assessing the efficacy of indirect methods.The current study aims to assess the performance of several indirect methods using LES data, which provide reliable velocity, turbulent stress, and bed shear stress data.Through these comparisons, this study aims to identify the most effective indirect method for complex flow fields.
The remainder of this paper is organized as follows.Section 3 provides a detailed description of the LES data and includes a comparative analysis of the LES-computed vertical profiles of the mean velocity and turbulent stress with ADV measurements for additional validation of the LES data.Section 4 presents comparisons of the bed shear stress calculated using the logarithmic law, RSS, and TKE methods, as well as their comparisons with the LES results.Finally, Section 5 concludes the paper with a summary of the major findings.

LES Model
The numerical method developed by Kang et al. (2011) is employed to conduct LES.This model solves the spatially filtered incompressible Navier-Stokes and continuity equations, closed with the dynamic Smagorinsky model, formulated in generalized curvilinear coordinates on a parallel supercomputer.The governing equations in the generalized curvilinear coordinates expressed using the Einstein summation convention are written as (i, j = 1,2,3) Water Resources Research 10.1029/2023WR036153 where J is the Jacobian of the geometric transformation given by J = |∂(ξ 1 , ξ 2, ξ 3 )/∂(x 1 ,x 2, x 3 )|, ξ i l = ∂ξ i /∂x l denote the transformation metrics, u i is the Cartesian velocity vector, U i = ξ i m /J) u m is the contravariant volume flux vector, g jk = ξ j l ξ k l is the contravariant metric tensor, p is the pressure, ρ is the density, μ is the dynamic viscosity, and τ ij is the sub-grid stress tensor arising from spatial filtering in the LES model.
The governing equations are solved using the parallel LES solver developed by Kang et al. (2011).This solver employs a matrix-free Newton-Krylov method for the momentum equation and an algebraic multigrid method for the Poisson equation.This approach enables fast convergence of the numerical solution, even on highly stretched grids often encountered in high Reynolds number flow modeling.Kang et al. (2021) conducted LESs for a laboratory experiment conducted by Jeon et al. (2018), which involved ADV measurements of velocity and turbulence statistics near and downstream of a spur dike, a type of hydraulic structure used for riverbank protection and river restoration.Kang et al. (2021) performed computations using approximately 57 million and 160 million grid points.The results of both computations were nearly indistinguishable, indicating that grid-converged solutions were obtained.Furthermore, these grid-converged results strongly agreed with the measured profiles of the mean velocity and Reynolds stresses.These high-resolution LES data, which employed a computational grid spacing small enough to resolve the viscous sublayer, were supported by validation results against experimental data, making them an ideal data set for evaluating the accuracy of indirect methods.

Description of the LES Data
In this study, the flow field data and bed shear stress distributions were obtained from LES data using 160 million grid points.Because this LES employed a minimum spatial resolution of 0.8 wall units, the wall shear stress could be calculated directly from the time-averaged flow field using Newton's law of viscosity, as presented in Equation 7, thus eliminating the requirement for ad hoc wall modeling.
In Equation 7, τ i and u i indicate the components of the wall shear stress and velocity in the i direction, respectively; and x, y, and z denote the streamwise, transverse, and vertical coordinates, respectively.The shear stress magnitude is calculated as Table 1 summarizes the flow and geometric parameters used for the LES and laboratory experiments (Jeon et al., 2018).The table indicates that the Reynolds number was moderately high, resulting in a fully turbulent flow.The Froude number was low; consequently, the effects of free-surface deformation and backwater were minimal, as reported by Jeon et al. (2018).
Figure 1 illustrates the computational domain.In the LES, a fully developed turbulent flow obtained from a precursor simulation entered the inlet (x/L = 6) and subsequently moved past a nonsubmerged spur dike (x = 0), which was attached to the right sidewall (y = 0) when viewed from the upstream side.Recently, Kang et al. (2023) discovered that the three-dimensionality of flow around spur dikes is governed by a wake parameter defined as S = c f L/H, where c f is the friction factor.They found that when the S value exceeded a critical value of approximately 0.04, the flow transitioned to being two-dimensional; below this value, it remained threedimensional.In the present LES, the wake parameter value was approximately 0.01.At this wake parameter value, the flow around a spur dike featured strong three-dimensional flow features, such as horseshoe and freesurface vortices, and recirculating flow regions in the wake and between the downstream face of the spur dike and the bed (Kang et al., 2021(Kang et al., , 2023)).Kang et al. (2021) already performed a comprehensively validated the LES results against the ADV measurements by comparing the transverse profiles of the computation and measurement at various cross-sections.Nonetheless, additional validations were presented to further corroborate the reliability of the LES data.In this section, additional comparisons are made between the LES computations and ADV measurements for the vertical profiles of the velocity and Reynolds stresses, as indirect methods estimate the bed shear stress based on these data.Before presenting the validation results, the variables and notation used in the remainder of this study are defined.u, v, and w refer to the instantaneous velocity components in the x, y, and z directions, respectively, and ρ implies the water density.The superscript ′ represents a temporally fluctuating component of a flow variable.τ LES refers to bed shear stress computed by LES, and u * = ̅̅̅̅̅̅̅̅̅̅̅̅ ̅ τ LES /ρ √ is the shear velocity.ν is the kinematic viscosity.U 0 indicates the mean bulk velocity.Additionally, k = (〈u′u′〉 + 〈v′v′〉 + 〈w′w′〉)/2 indicates TKE.
Figure 2 presents a comparison between the LES and measurements at locations of x/L = 0.83 and 1.67.This comparison involves the vertical profiles of 〈u〉, 〈v〉, ρ〈u′u′〉, ρ〈w′w′〉, and ρ〈u′w′〉, all of which were made dimensionless by U 0 .The computed 〈u〉 and 〈v〉 values are in strong agreement with the ADV measurements.
Although not presented here owing to space constraints, the computed 〈w〉 values also showed a comparable level of agreement.In summary, the LES reproduced the measured velocity field with good accuracy.
Additionally, the computed ρ〈w′w′〉 and ρ〈u′w′〉 Reynolds stresses exhibited good agreement with the ADV measurements.However, there are some discrepancies in the ρ〈u′u′〉 stress.This can be attributed to the wellknown overestimation of the horizontal velocity fluctuations in ADV.For instance, Voulgaris and Trowgridge (1998) showed that horizontal velocity fluctuations-〈u′u′〉 and 〈v′v′〉 in our case-as measured by ADV, were prone to high-frequency noise, resulting in an overestimation of these values.Nonetheless, the vertical velocity fluctuations, mean velocities, and Reynolds shear stresses measured using the ADV were in excellent agreement with those obtained using laser Doppler velocimetry.Similar observations were reported by Khorsandi et al. (2012).Through jet experiments, the authors demonstrated that the ADV overestimated the horizontal velocity fluctuations when compared to measurements taken using hot-film anemometry, whereas the vertical velocity fluctuations exhibited satisfactory agreement.They attributed this overestimation to the Doppler noise and spikes.The two studies strongly suggest that the discrepancy in the ρ〈u′u′〉 stress shown in Figure 2 could be linked to the inherent overestimation of horizontal velocity fluctuations in the ADV measurements.This hypothesis is further supported by the strong agreement observed in all other flow quantities such as 〈u〉, 〈v〉, 〈w〉, ρ 〈w′w′〉, and ρ〈u′w′〉, with the discrepancy being evident only in the ρ〈u′u′〉 profile.
In summary, the present validation results showed a favorable agreement between the LES and ADV measurements.This result, in conjunction with the additional validation and grid sensitivity tests previously conducted by Kang et al. (2021), confirmed the accuracy and grid convergence of the LES results.In Section 4, the LES-computed flow field data are employed to assess the performance of various indirect methods.

Results
The bed shear stresses calculated from the LES results using Equation 7, henceforth referred to as τ LES , served as a reference value for evaluating the accuracy of the indirect methods examined in this study.In the LES data, there were approximately 7.3 × 10 5 grid points at the bottom plane.

10.1029/2023WR036153
As illustrated in Figure 1, the flow domain was divided into three regions: Regions 1, 2, and 3. Region 1, located upstream of the obstacle (x ≤ 2L), represents a region characterized by a predominantly unidirectional mean flow.Region 2, located near the obstacle ( 2L < x ≤ 3L), encompasses the near-wake region where energetic vortex shedding, flow separation, and intense turbulent mixing are observed.Last, Region 3 comprises the remaining wake region (3L < x ≤ 23.33L), where streamlines separating from the spur dike's tip reattach to the sidewall, forming a large recirculation zone elongated in the streamwise direction.In Region 3, the turbulence levels were considerably lower than those in Region 2, and the reduced velocity in the wake region gradually recovered as the flow moved away from the obstacle.
In the subsequent subsections, the results of the logarithmic law, RSS, and TKE methods are presented and compared with reference shear stress values.

Logarithmic Law
The logarithmic law assumes that when the velocity and distance from the wall are both made dimensionless using the shear velocity, the plot of the dimensionless velocity against the dimensionless distance exhibits a linear trend on a semi-logarithmic chart.This implies that if this law is valid, a distinct linear section should exist in the plot.Figure 3 shows the magnitudes of the horizontal velocities (〈u〉 and 〈v〉) against z + .The dashed lines in the figure represent Equation 2, with constant values of κ = 0.40 and B = 5.
Figure 3a shows the distinct linear sections within z + < 10 3 in Region 1.A decrease in velocity that is noticeable slightly below the free surface is associated with the velocity dip phenomenon caused by sidewall effects (Yang et al., 2004).The existence of this logarithmic region implies that in Region 1, the bed shear stress can be accurately estimated using the logarithmic law with good accuracy.A notable finding is that the logarithmic region diminished in size as the section approached an obstacle (x = 0).Figure 3b, presenting results for Regions 2 and 3, exhibits markedly different patterns.In sections where 0 < x ≤ 5L, the velocity plots display pronounced irregularities, with the logarithmic region being absent in most sections.This absence of the logarithmic region is possibly owing to the complex and three-dimensional velocity patterns around the spur dike, influenced by recirculating flows and the horseshoe vortex system.In further downstream sections where x ≥ 10L, the logarithmic section starts to reappear in the velocity plot, indicating that the flow gradually reverts to a one-dimensional state in the downstream direction.
These findings indicate that the logarithmic law is effective only in Region 1, where the flow is quasi-onedimensional with nearly straight streamlines.Furthermore, the results showed that the logarithmic region may be completely absent when strong three-dimensional flow and wake effects exist.These results indicate that the logarithmic law method may be unsuitable for estimating bed shear stress in areas influenced by obstacles.This limitation arises from the fact that this law is applicable only to simple boundary layer flows.To enhance bed shear stress estimation in regions where the logarithmic law cannot be applied, indirect methods utilizing RSS or TKE data have been proposed.

RSS Method
The underlying assumption behind the RSS method is that an RSS value close to the bed is a good approximation of the local bed shear stress (Graf & Song, 1995;Song & Chiew, 2001;Song et al., 1994).Various techniques for choosing the near-bed value have been used in past research, such as the single-point, linear extrapolation, and maximum stress methods.The bed shear stress values estimated using the single-point, linear extrapolation, and maximum RSS methods are denoted as τ RSS(single) , τ RSS(ext) , and τ RSS(max) , respectively.τ RSS(ext) was further split into τ RSS(ext,A) and τ RSS (ext,B) , each employing different extrapolation points.Further details are provided as follows.
First, the single-point method employs one point from an RSS profile located close to the bed to estimate the local shear stress.This method has been predominantly used in oceanographic studies (Heathershaw, 1979;Heathershaw & Simpton, 1978;Huthnance et al., 2002;Kim et al., 2000), with the chosen locations differing between studies.Heathershaw (1979) and Heathershaw and Simpton (1978) selected a point 1.0 m above the bed.Conversely, Huthnanceet al. (2002) and Kim et al. (2000) used points at 0.30 and 0.14 m above the bed, respectively.However, Biron et al. (2004)  approximately 0.1 times the flow depth above the bed.In this study, τ RSS(single) was calculated using various options, including those suggested by Biron et al. (2004).For Region 1, heights of z = 0.02 and 0.10H were chosen, whereas for regions 2 and 3, heights of z = 0.02, 0.06, and 0.10H were used.
The linear extrapolation method determines the bed shear stress by extending an RSS profile to the bed and then identifying the RSS value at which it intersects (Biron et al., 2004;Graf & Song, 1995;Nikora & Goring, 2000;Yang et al., 2004).There are different approaches to this extrapolation depending on the shape of the RSS profiles.For example, not all RSS profiles have distinct linear sections, particularly when shear layers are present.Lichtneger et al. (2020) highlighted that the method of extrapolation can often overestimate the bed shear stress in nonuniform flows because the maximum RSS may be located considerably away from the bed.Typically, the rise in turbulent stress, including the RSS, results from a combination of turbulence from the shear layers and bottom friction (Shiono & Muto, 1998).The turbulence intensity from the shear layers could be orders of magnitude greater than that from bottom friction, as demonstrated in our earlier LESs and experiments on flows past obstacles in open channels (Jeon et al., 2018;Kang et al., 2016Kang et al., , 2021Kang et al., , 2023)).The flow field considered in the present work features shear layers induced by vortex shedding emanating from the tip of the spur dike and flow separation around and downstream of the obstacle.In these cases, determining ways to extrapolate the RSS profile to the bed is ambiguous.Two extrapolation techniques are evaluated in this section.τ RSS(ext,A) was calculated by performing a linear regression on a linear segment of the RSS profile, typically within 0.2H < z < 0.8H.Conversely, τ RSS(ext, B) was determined by forming a straight line using two points within z ≤ 0.10H.τ RSS(ext,A) was calculated for Region 1, where distinct linear segments are readily identifiable in the RSS profile.For Regions 2 and 3, τ RSS(ext,B) was used Regions 2 and 3 because defining the linear segments can be difficult.
Finally, the maximum shear stress method considers the maximum value within the RSS profile to be the bed shear stress.As mentioned previously, the presence of shear layers can shift the RSS peak position upward.When this occurs, the method can significantly overestimate the bed shear stress.For instance, Lichtneger et al. (2020) reported that in these cases, a linear extrapolation method can yield an excessively large shear stress, making the maximum RSS method a better choice than the extrapolation method.
A crucial consideration in determining bed shear stress using RSS methods that is often overlooked is the inclusion of the transverse RSS component.In most previous studies, only the primary RSS component, 〈 ρu′w′〉 in our case, was considered.However, in our case, the transverse component of the bed shear stress can be significant in the regions around the spur dike, where the flow direction shifts away from the primary flow direction, or inside the downstream recirculation region, when the velocity magnitude is considerably less than the bulk velocity.To address these effects, the transverse RSS component was considered when calculating bed shear stress values in this study.This bed shear stress magnitude is defined as Although the inclusion of the transverse component might have had a minimal effect in Region 1, it was significant in Regions 2 and 3.

Region 1
Figure 4 displays the vertical profiles of dimensionless RSS components in the x and z directions, represented as ⟨ u′w′⟩/ U 2 0 and ⟨ v′w′⟩/ U 2 0 , respectively, for Region 1.The 〈 u′w′〉 profile exhibited a local minimum at the velocity dip location (z ˜0.9H).From this point, its value rises linearly to the peak at z ˜0.1H.Following the peak, the value diminishes to zero as it approaches the bed.This pattern is consistent with the findings of other experimental studies (Biron et al., 2004;Graf & Song, 1995;Nikora & Goring, 2000;Yang et al., 2004).Conversely, the magnitude of 〈 v′w′〉 is significantly smaller than that of 〈 u′w′〉 across all locations.This can be attributed to the flow direction in Region 1 being predominantly aligned in the x-direction.In this region, the contribution of the transverse RSS component to the magnitude of the bed shear stress is likely to be minimal.Biron et al. (2004) employed a single-point method using the points at z = 0.10H, which seems to be an appropriate selection given that the peak's location for 〈 u′w′〉 varied between 0.04 and 0.10H, as shown in Figure 4. To assess the effect of the sampling location on the results, τ RSS(single) was computed using two different heights, z = 0.02 and 0.10H.Moreover, τ RSS(ext,A) was determined by fitting the RSS profile between z/H = 0.2 and 0.7 to a linear equation and subsequently extending it to the bed.τ RSS(max) was estimated using the maximum RSS value in the profile.

10.1029/2023WR036153
Figure 5 presents comparisons of τ RSS(ext,A) , τ RSS(single) , and τ RSS(max) with τ LES along with their respective root mean square errors (RMSEs) relative to τ LES .The sequence of smallest RMSEs is τ RSS(max) , τ RSS(single) , and τ RSS (ext, A) .However, differences between the results were minimal.In Figure 5b, the single-point method using z = 0.1H exhibited a marginally smaller error than that using z = 0.02H.These results indicate that for Region 1, where the flow direction is aligned with the primary flow direction, (a) the extrapolation method is slightly less accurate than the other two methods, (b) the single-point method's accuracy is not strongly influenced by the sampling location, and the one with z = 0.1H employed by Biron et al. ( 2004) is a reasonable choice, (c) both the single-point and Water Resources Research 10.1029/2023WR036153 maximum RSS methods yield similar results, and (d) despite some discrepancies, the overall accuracy of all three RSS methods for Region 1 is encouraging.

Regions 2 and 3
The flow patterns in Regions 2 and 3 are significantly more complex than those in Region 1 owing to the interactions between the flow and the spur dike.This results in recirculating flow regions, strong secondary flows, energetic vortex shedding, and increased turbulence intensity (Jeon et al., 2018;Kang et al., 2021Kang et al., , 2023)).Typically, the velocity profiles in these regions do not exhibit distinct logarithmic sections, which makes the application of the logarithmic law challenging.
Figure 6 shows the computed RSS profiles in sections between 2L ≤ x ≤ 20L.This figure presents markedly different RSS profiles, compared to those for Region 1 in Figure 4. First, the magnitude of 〈 v′w′〉 is no longer negligible and comparable to that of 〈 u′w′〉.For instance, at the position x = 3L and y = 1.33L, which lies within the extent of the shear layer formed by vortex shedding from the tip of the spur dike, the peak value of 〈 v′w′〉 is approximately 60% of 〈 u′w′〉.This suggests that the inclusion of this transverse RSS component may be crucial for estimating the bed shear stress in these regions.Second, the 〈 u′w′〉 profiles no longer exhibit a distinct linear section, and the vertical positions of the peak RSS values vary considerably between locations.The absence of linear sections together with the irregular positioning of peak RSS values can present challenges in applying linear extrapolation and maximum stress methods.
Figure 7 shows the distributions of the peak RSS vertical coordinates, denoted as z max , along with their probability distributions throughout the flow domain.Figure 7a shows that in the majority of Region 1, the z max /H values are less than 0.1H.However, in areas close to the spur dike and the downstream recirculating flow region, z max /H exceeds 0.5.Beyond the reattachment point (x ≈ 12L), z max /H decreases in magnitude, and its distribution becomes more uniform in the lateral direction.These results indicate that z max /H can fluctuate considerably depending on the local flow characteristics.In Figure 7b, the peak RSS values occurring within z < 0.1H account for approximately 80%, 50%, and 15% of the data points in Regions 1, 2, and 3, respectively.Arrows in this figure highlight the peaks of the probability distributions across Regions 1-3, signifying the characteristic z max values for each area.Interestingly, Regions 1 and 2 share similar peak locations (z max /H = 0.02-0.04),diverginf sharply from Regions 3, which shows a markedly higher value (0.24).These differences are likely due to the different flow patterns between Regions 3 and the others.Specifically, Region 3, situated in the far wake of the spur dike, exhibits complex flow patterns due to the interactions of spatially evolving multiple turbulent shear layers.Such pronounced variability in z max values, which are influenced by local flow patterns, could present a difficulty in the determination of optimal sampling locations for RSS methods.These methods conventionally assume a strong correlation between peak RSS values and bed shear stress.Subsequently, the influence of the sampling location on the estimated shear stress is examined for both single-point and linear extrapolation techniques.

10.1029/2023WR036153
Figure 8 presents the bed shear-stress distribution estimated using linear extrapolation, single-point, and maximum stress methods juxtaposed with τ LES .For the single-point approach, the results using z/H values of 0.02, 0.06, and 0.10 were plotted.For the extrapolation approach, the results with combinations of z/H = 0.02 and 0.04, z/H = 0.06 and 0.08, and z/H = 0.10 and 0.12 were displayed.As shown in Figure 8b, the linear extrapolation method exhibited satisfactory agreement with the LES results.The τ RSS(0.02H-0.04H ) showed excellent agreement with τ LES .However, τ RSS(0.08H-0.10H)somewhat overestimated the values, compared to the other two cases.This suggests that the accuracy improves when the sampling point, particularly the lower point, is closer to the bed.In

10.1029/2023WR036153
Figure 8c, the τ RSS(single) distributions are in reasonable agreement with τ LES in Region 1.However, in Regions 2 and 3, notable differences emerge, compared to the LES result, particularly for z = 0.06 and 0.08H.For these two cases, the τ RSS(single) distribution features a thin layer of low shear stress surrounded by zones of increased stress.This layer originated from the tip of the obstacle and extended toward the y = 0 sidewall.However, this distribution is not evident in the LES results, as shown in Figure 8a.Among the three sampling positions, z = 0.02H yielded the most qualitatively similar results to the LES data.Figure 8d shows the result of the maximum RSS method.A clear discrepancy was observed when juxtaposed with the τ LES distribution.This method overestimates the extent of the high shear stress zones around the tip of the spur dike.This observation was somewhat expected, given that the shear layers around the spur dike can increase the RSS values at locations distant from the bottom.
In summary, the linear extrapolation method yielded the most promising results among the three indirect methods, whereas the maximum stress method was the least accurate.It is worth noting that all RSS methods struggle to predict the shear stress patterns near the tip of the spur dike.In particular, they failed to capture the presence of a thin crescent-shaped layer with increased shear stress in the region circled in Figure 8a.According to Kang et al. (2021Kang et al. ( , 2023)), this layer is the footprint of a horseshoe vortex system moving around a spur dike.This result suggests that indirect methods may be unable to capture the complex physics associated with the horseshoe vortex system around obstacles.The above results clearly demonstrate that the two-point extrapolation method outperforms the other two RSS methods.This leads to the inquiry: Which combination of sampling points yields the most accurate result, and what is the reason behind this outcome?To explore this, the RMSEs of extrapolation methods with various combinations of sampling positions are examined in Table 2. Additionally, for comparison purposes, the RMSEs for both single-point and maximum stress methods are listed.
Table 2 indicates that the total RMSE, representing the RMSE for the entire domain, as well as the RMSE in each region, is the smallest in Cases E12 and E13, which use z = 0.02H as the lower sampling point.Other cases that use the sample lower sampling point but with different upper points, namely Cases E11, E14, and E15, also show similarly favorable results.Additionally, it can be observed that other RSS methods utilizing extrapolation, Cases E21-E55, exhibit increased RMSEs as the vertical coordinates of the lower sampling points increase.Another notable finding is that the accuracy of the single-point method improves as the sampling position approaches z = 0.02H.However, in comparison to extrapolation methods employing the same lower sampling point as the single-point method, the extrapolation approach generally outperforms the single-point method.Lastly, the maximum stress method exhibits the least favorable performance among all the methods evaluated.
Several key findings emerge from the foregoing results.First, the extrapolation method demonstrates superior efficacy compared to both the single-point and maximum stress methods.Second, within the extrapolation method, selecting the lower sampling point proves to be more critical to achieving accuracy than the choice of the upper point, with z = 0.02H as the lower point yielding the most favorable outcomes.Lastly, extrapolation methods employing this particular lower point consistently yield some of the lowest RMSEs across all sub-  regions, Regions 1-3, despite the significant variability observed in peak RSS positions, as depicted in Figure 7.This suggests that the optimal location of the lower sampling point may not be sensitive to local flow characteristics, such as peak RSS positions.For instance, Region 3, with a characteristic z max /H value of approximately 0.12 (Figure 7), shows decreasing accuracy as the lower sampling location approaches this height.
Figure 9 shows a quantitative comparison of the bed shear stress predictions between the RSS method and the LES data.In Figure 9a, the results of the extrapolation methods with the lower points at z = 0.02, 0.06, and 0.10H are presented.It is important to note that the spacing between the upper and lower points is 0.02H in Cases E11, E31, E51, 0.06H in Cases E13, E33, and E53, and 0.10H in Cases E15, E35, and E55.This figure reveals a distinct trend where both RMSE and R 2 values increase, indicating a deterioration in accuracy, as the lower sampling point is positioned further from the bed.Additionally, the effect of the sampling point spacing on the errors seems to be insignificant.Figure 9b displays the results of single-point methods.The comparisons between Cases E11 and S1, E13 and S3, and E15 and S5 demonstrate that, when applied at identical elevations, the extrapolation method consistently outperforms the single-point method in terms of both RMSE and R 2 values.In Figure 9c, the  9d compares the RMSEs in the sub-regions (Regions 1, 2, and 3) for Cases E12, S1, and M. In Region 1, all cases yield similar outcomes with low RMSEs, suggesting that any RSS method can be effectively applied to straight-flow regions.However, their performances diverge significantly in Regions 2 and 3.In these region, the flow patterns are complex and highly three-dimensional, characterized by significant flow separation from the spur dike's tip and the horseshoe vortex systems in front of it.This indicates that the extrapolation method is well-suited to regions with complex flow patterns.Another important observation is that the shear-stress predictions using only the streamwise RSS component have a higher RMSE than the total stress accounting for both streamwise and transverse RSSs in the extrapolation and single-point methods.This indicates that incorporating the transverse RSS can improve prediction accuracy, particularly in complex flow scenarios.For instance, in Cases E12 and S1, the inclusion of transverse shear stress in the shear stress calculations resulted in accuracy improvements of 8% and 3% within Region 2, as measured by RMSE values, respectively.These improvements were more pronounced in Region 2*, a 2L × L rectangular zone adjacent to the tip of the spur dike.The flow in this region is significantly influenced by the deflection of upstream flow toward the channel center, resulting in strongly curved streamlines.
In this specific region, accuracy improvements of 13% and 3% were observed for Cases E12 and S1, respectively.Such findings underscore the benefits of including transverse shear stress in enhancing the accuracy of the extrapolation methods in areas with strong secondary flows.
The TKE method is primarily based on observations of coastal boundary layer flows, suggesting a linear relationship between bed shear stress and TKE close to the bed (Stapleton & Huntley, 1995).Soulsby (1983) indicated that the proportionality constant C 1 in Equation 4 typically lies between 0.18 and 0.21.Other researchers (Hopkinson & Wynn-Thompson, 2012;Stapleton & Huntley, 1995) proposed a constant value of 0.19.This value has been extensively used in oceanographic contexts (Pope et al., 2006;Stapleton & Huntley, 1995;Thompson et al., 2003).However, some researchers noted that this value, which originates from tidal flow measurements, may not be suitable for riverine flows (Biron et al., 2004;Hopkinson & Wynn-Thompson, 2012;Noss et al., 2010).Therefore, some studies have attempted to evaluate TKE methods for fluvial applications.For instance, using streamflow data, Noss et al. (2010) demonstrated that TKE and RSS data do not exhibit a linear relationship with a unique C 1 value.Hopkinson and Wynn-Thompson (2012) conducted a laboratory experiment to determine the C 1 value, noting that earlier studies derived C 1 values by using RSS values rather than directly measuring bed shear stress.They performed a direct measurement of bed shear stress, and from these data, C 1 values ranging from 0.11 to 0.53 were obtained, suggesting potential inaccuracies when TKE methods were used for stream flows.Some studies have questioned the applicability of TKE methods to fluvial applications.
In this section, the efficacy of the TKE method was evaluated using the LES-computed velocity field and bed shear stress data.In Figure 10, k and τ LES are plotted for points at z = 0.02 and 0.1H.From the plots, C 1 values can be estimated from the slope of a best-fit line for the flow data using the relation C 1 = τ LES /(ρk).The solid black line represents the line with constant C 1 = 0.19 typically used in oceanographic research, and the dashed lines represent the best-fit lines for the data sets from the two sampling heights.
At x/L = 5 and 4, a linear relationship between k and τ LES is visible.For a sampling height of 0.02H, C 1 values are within the 0.21 to 0.22 range, and they fall within the 0.18-0.21range proposed by Soulsby (1983).For a sampling height of 0.10H, C 1 values increase to the 0.28-0.30range.These results show the existence of a linear correlation between the TKE and bed shear stress in the upstream straight-flow region.In this region, the C 1 values at the sampling location of z = 0.02H were similar to those observed in the tidal boundary-layer flows.This implies that tidal boundary-layer flows exhibit flow characteristics similar to those of a boundary layer flow in a straight channel.within the recirculation zone.The data points were most dispersed at x = 3L, with the dispersion diminishing in the x direction.Downstream of x = 12L, which is the reattachment point of the recirculation zone, the relationship between k and τ LES starts to revert to a linear pattern.Interestingly, at x = 22L, the scatterplots show a linear pattern defined as τ LES = C 1 k + α, with α being a nonzero constant.By comparing the results at x = 18 and 22L, it becomes evident that the value of α approaches zero, and the slope tends toward 0.19 as the flow recovery continues.
The above results show that a linear relationship between the TKE and bed shear stress, as suggested by various researchers, is valid only under limited conditions.This linear relationship was most evident in the areas where the Figure 11 shows the distributions of τ LES , k, and τ LES /(ρk) in the z = 0.02H plane.τ LES and k show similar distributions in x < 3L, and in this region, τ LES /(ρk), corresponding to a local value of C 1 , exhibits an almost constant value.In other areas, the τ LES and k distributions did not show a good correlation.In addition, the τ LES / (ρk) values displayed significant variations.They tended to be low (<0.05) within the recirculation region but exceeded 0.4 near the obstacle.This result demonstrates that a universally applicable C 1 value does not exist in complex open-channel flows.

Summary and Conclusions
This study assessed the performance of indirect methods for predicting bed shear stress.Three indirect methods were tested for an open-channel flow past a spur dike: the logarithmic law, RSS, and TKE methods.Highresolution velocity and turbulent stress data from our previous LES were employed for application and comparison of these indirect methods.The main conclusions of this study are as follows: The velocity data exhibited a distinct logarithmic section upstream of the spur dike.However, this section was less evident in other areas, where the flow direction was not well aligned with the primary direction.This suggests that the logarithmic law is best when applied to flows characterized by nearly straight streamlines and a limited transverse velocity.
Second, the single-point, two-point extrapolation, and maximum stress RSS methods showed promising results in areas upstream of the spur dike where the flow was not affected by the spur dike.However, in areas influenced by the spur dike, the results of the three methods varied considerably.Among the three methods, the two-point extrapolation method yielded the most accurate results, particularly when the lower sampling points were closer to the bed.The single-point method also performed reasonably when the sampling point was close to the bed.Conversely, the maximum stress method was the least accurate because it interpreted the peak RSS of the turbulent shear layers as the bed shear stress.In addition, the inclusion of the transverse RSS component significantly improved the accuracy of the calculated shear stress, particularly around the spur dike, where the flow patterns were the most complex.The above observations lead to a practical recommendation: the accuracy of the single-point method, commonly employed at a sampling height of 0.1H, can be improved by introducing an additional sampling point either below or above this height and employing linear extrapolation, without incurring significant extra effort.
Third, the TKE method proved effective in areas upstream of the spur dike where a constant value of C 1 existed.In other areas, there was no linear relationship between TKE and bed shear stress.The C 1 value, inferred from the distributions of bed shear stress and TKE, showed significant variability around and downstream of the spur dike.It ranged from 0 to 0.5.Therefore, it is not recommended to use TKE methods for flows involving obstacles.In addition, the C 1 value upstream of the spur dike showed a clear dependence on sampling height.This suggests that employing the TKE method, even for simple open-channel flows, may pose challenges.
In conclusion, the RSS method, using a two-point extrapolation technique with a lower point located close to the bed, has been proven to be effective in all flow regions and outperformed the TKE and logarithmic law methods.
A major limitation of this study is that these conclusions might be specific to the flow conditions analyzed here, which are characterized by low Froude numbers, low Reynolds numbers, and smooth-wall conditions.Consequently, additional research is essential to conduct further research to establish an extrapolation method capable of yielding optimal results under a wide range of flow conditions.
To the best of our knowledge, this study is the first to assess the accuracy of bed shear stress predictions by indirect methods in a large flow domain involving an obstacle, by comparing them with direct shear stress measurements.Subsequent investigations could explore flows under various flow and geometric conditions, as well as examine the potential suitability of these techniques for natural streams.

Figure 1 .
Figure 1.Computational domain for the large-eddy simulation.The x, y, and z coordinates indicate the streamwise, transverse, and vertical directions, respectively.

Figure 5 .
Figure 5.Comparison of the bed shear stress estimated using Reynolds shear stress methods for Region 1 with the large-eddy simulation data.Points spaced 0.11H (0.1 m) apart in both x and y directions, resulting in 89 points in the streamwise direction and eight points in the transverse directions, are used for the analysis.

Figure 7 .
Figure 7. (a) Vertical coordinates (z max /H) of the peak Reynolds shear stress and (b) their probability distributions.
At x = 3L, the linear relationship between k and τ LES begins to break down.At x = L, it disappears completely.C-shaped curves emerge in the downstream sections.The lower segment of the curve corresponds to the area

Figure 10 .
Figure 10.Scatterplot showing a correlation between τ LES and k at z = 0.02 and 0.10H.S denotes a slope.

Table 1
Flow and Geometric Parameters Used in the Large-Eddy Simulation

Table 2
Root Mean Square Errors for Reynolds Shear Stress Methods Comparison of bed shear stresses estimated using Reynolds shear stress methods with the large-eddy simulation data.In panel (d), the total and streamwise (xcomponent) shear stresses are compared across different sub-regions.R 2 is the square of linear correlation.Points spaced 0.11H (0.1 m) apart in both x and y directions, resulting in 89 points in the streamwise direction and eight points in the transverse directions, are used for the analysis.Region 2* is a rectangular zone bounded by 1 ≤ x/L ≤ 1 and 1 ≤ y/L ≤ 2. stress method yields the least accurate results, as discussed earlier.Lastly, Figure maximum