Bed Shear Stress and Near‐Bed Flow Through Sparse Arrays of Rigid Emergent Vegetation

Vegetation is an essential component of natural rivers and has significant effects on flow and morphodynamic processes. Although progress has been made in characterizing flow resistance in vegetated flows, the impact of vegetation on bed shear stress, a key driver of sediment transport, still needs better characterization and understanding. This research, explores bed shear stress and near‐bed flow characteristics in sparse arrays of rigid emergent cylinders mimicking vegetation over a rough bed. For this purpose, a novel adaptation of a shear plate was used to measure bed shear stress at the canopy scale. These measurements were analyzed in relation to spatially averaged near‐bed flow quantities for different array densities. The results show that, for a constant water depth, the investigated cylinder canopy enhances the ratio between bed shear stress and bulk flow velocity (i.e., Darcy‐Weisbach bed friction factor) compared to unobstructed open‐channel flows, and that this ratio increases with array density. Moreover, higher near‐bed velocities were observed for higher array densities. On the other hand, no influence of the cylinder array on near‐bed values of turbulent kinetic energy and turbulent stresses was observed. Finally, it is shown that the thickness of the near‐bed layer is a suitable parameter to scale the ratio between bed shear stress and bulk flow velocity in vegetated rough bed flows.


Introduction
Riparian and in-stream aquatic vegetation are essential components of fluvial systems, and their interplay with flow and sediment transport processes is considered beneficial for the health status of fluvial ecosystems (e.g., Bywater-Reyes et al., 2022).The interaction of vegetation with flow and sediment processes alters the flow field and sediment transport rates at different spatial scales.From an ecological perspective, these changes improve the environmental conditions for sustaining aquatic life by enhancing the formation of heterogeneous habitat patches and improving water quality by trapping harmful substances (Rowiński et al., 2018).On the other hand, vegetation can increase flood risk by raising water levels due to its contribution to flow resistance and it can enhance the risk of failure of hydraulic structures due to significant changes in erosion and deposition rates.
Current bed load transport models developed for vegetated flows require accurate estimates of bed shear stress (e. g., Armanini & Cavedon, 2019;Wang et al., 2023) or near-bed turbulent kinetic energy (e.g., Zhao & Nepf, 2021) for their application.Although significant efforts have been made to improve the understanding and modeling of the impact of vegetation on flow resistance (e.g., Box et al., 2022;D'Ippolito et al., 2023) and turbulence production (e.g., Ricardo et al., 2014;Xu & Nepf, 2020), there exist still knowledge gaps related to the estimation of bed shear and near-bed flow quantities within vegetated areas (D'Ippolito et al., 2023).Therefore, proper characterization of the influence of vegetation on bed shear stress and the near-bed flow structure is a critical step to improve morphodynamic modeling and elucidate the particular mechanisms governing sediment transport in vegetated flows (Wu et al., 2021).
Although several methods are available for the determination of bed shear stress in non-vegetated open-channel flows such as the depth-slope -, near-bed Reynolds stress-, and turbulent kinetic energy method (Biron et al., 2004;Rowiński et al., 2005), they cannot be directly applied to vegetated flows due to the interaction of the vegetation with the flow (Le Bouteiller & Venditti, 2015;Nepf, 2012).Therefore, different methods have been specifically developed for vegetated flows (Conde-Frias et al., 2023;Yang et al., 2015) considering mainly smooth beds.In the case of vegetated rough bed flows, it is not clear how the bed shear stress relates to the nearbed flow structure.Recent measurements by Penna et al. (2020) carried out with emergent rigid vegetation and rough beds revealed that the interplay between vegetation and bed roughness significantly impacts near-bed streamwise velocities and turbulent kinetic energy.This shows the need for a better characterization of the combined impact of vegetation and bed roughness on the near-bed flow structure and bed shear stress.This paper investigates the impact of vegetation array density on bed shear stress and its relevance for near-bed flow quantities over a rough bed through the combination of bed shear stress measurements using a novel shear plate and detailed flow measurements in the near-bed flow region.The paper is structured as follows: Section 2 presents a review of the impact of rigid emergent vegetation on hydraulic resistance, near-bed flow structure, and its linkages with bed shear stress.Section 3 describes the experimental setup and procedure.Section 4 presents and discuss the results and Section 5 highlights the main findings of this paper.

Hydraulic Classification of Vegetation
Vegetation in riverine environments is diverse in terms of plant characteristics and array configurations.From a hydraulic engineering perspective, vegetation is commonly classified dependent on its degree of submergence, flexibility, and array density (DWA, 2020).The degree of submergence refers to the ratio between flow depth h and vegetation height h v , that is, vegetation is classified as emergent and submerged if h/h v ≤ 1 and h/h v > 1, respectively.Flexible vegetation adapts its height and shape in response to the flow conditions resulting in a variable frontal area exposed to the flow A v , whereas for rigid vegetation A v varies only with water depth (A v = d v h, where d v is the stem diameter for the case of cylindrical elements).Array density is commonly characterized by the stem areal density m = 1/Δs 2 (stems/m 2 ), where Δs = distance between stems, or as solid volume fraction occupied by the vegetation, ϕ v = md v 2 (π/4) (-).For the case of rigid emergent vegetation, the array density can be classified according to the diameter spacing ratio d v /Δs into sparse (d v / Δs < 0.56), and dense (d v /Δs ≥ 0.56) based on the characteristic scale governing turbulence production (Tanino & Nepf, 2008).The present research focuses on sparse arrays of rigid emergent vegetation.

Impact of Vegetation on Flow Resistance
In vegetated flows, flow resistance is governed by the drag force F d exerted by the flow on each plant and by the bed shear stress τ b (Figure 1).For arrays of rigid emergent vegetation, the spatially averaged drag force acting on the canopy can be parametrized following the "rigid cylinder analogy" (e.g., Aberle & Järvelä, 2013), as where ρ = water density, C d = mean bulk drag coefficient of the array, 〈h〉 = spatially averaged flow depth (angled brackets denote spatial averaging), U = bulk flow velocity, with U = Q/(B〈h〉), Q = flow discharge, and B = channel width.
The drag coefficient of an isolated cylindrical object C d depends on the stem Reynolds number Re s = Ud v /ν, where ν = kinematic viscosity, but is almost constant for 8 × 10 2 < Re s < 2 × 10 5 (e.g., Schlichting & Gersten, 2017).However, in the case of arrays, the bulk drag coefficient C d differs from the drag coefficient of an isolated element (Lindner, 1982) due to the heterogeneous flow field associated with cylinder arrangement and array density.Although different models have been developed to assess the impact of array density on C d (e.g., Etminan et al., 2017;Lindner, 1982;Petryk & Bosmajian, 1975;Sonnenwald et al., 2019), there exist still open questions regarding the effect of rigid emergent canopy characteristics on C d such as array density and cylinder arrangement (Aberle & Järvelä, 2013).Therefore, direct force measurements have been used for the determination of 〈F d 〉 and hence C d (e.g., Kothyari et al., 2009;Schoneboom et al., 2011;Tinoco & Cowen, 2013).

Impact of Vegetation on the Flow Field
The vertical flow structure over rough beds with rigid emergent vegetation can be subdivided into three layers as shown in Figure 1 (Penna et al., 2020;Ricardo et al., 2016).In the upper layer, the flow is mostly affected by the wakes induced by local flow disturbances at the water surface (z ws ) due to the emergent vegetation elements.This region is sufficiently far away from the bed so that it has no impact on the near-bed flow structure.In the intermediate layer (IL), flow resistance and turbulence production are governed by the stem wakes.For vertically Water Resources Research 10.1029/2023WR035879 uniform arrays, the flow velocity and turbulence parameters are constant over depth and the streamwise velocity is almost equal to the bulk flow velocity (Conde-Frias et al., 2023;Ricardo et al., 2016).In the near-bed layer (NBL), also called the bottom boundary layer (Conde-Frias et al., 2023), the flow field is highly threedimensional (3D) and is controlled by the interplay between vegetation and the bed roughness (Penna et al., 2020).Although the flow characteristics in the NBL govern bed shear stress, few studies have been carried out focusing on this layer (Conde-Frias et al., 2023;Yang et al., 2015Yang et al., , 2016)), especially for rough bed flows (Penna et al., 2020).
The analytical description of the spatially heterogeneous near-bed flow field in vegetated conditions can be achieved by using the double-averaging methodology (DAM; Nikora et al., 2007), that is, by averaging flow properties in time and space using thin slabs oriented parallel to the bed.The spatial and temporal decomposition of the flow quantities results in temporal and spatial fluctuation components θ′ = θ θ and θ = θ ⟨θ⟩, respectively, where θ is the instantaneous value of the flow quantity, θ the temporally averaged value and ⟨θ⟩ the double-averaged (in time and space) quantity.With these definitions, the temporally and spatially averaged momentum equations (double-averaged Navier-Stokes equations [DANS]) can be derived for steady quasiuniform flows using the Einstein notation (Ferreira et al., 2009;Nikora et al., 2004) according to where g i = ith component of the gravity acceleration; ϕ w = fluid volume fraction within the averaging slab; u i = ith component of the velocity vector; f v = sum of viscous and pressure drag per unit fluid volume due to vegetation; f b = sum of viscous and pressure drag per unit fluid volume due to bed roughness.The first and second terms on the right-hand side stand for the viscous and turbulent (Reynolds) stresses respectively.The third term represents the so-called form-induced stresses.
Using DAM, analytical expressions for the flow field in the IL can be derived.Assuming that flow resistance is governed by the drag forces and that vertical gradients of turbulent and form-induced stresses are negligible, the double-averaged streamwise velocity in the IL 〈u〉 IL can be formulated according to Nikora et al. (2004) 〈u〉 where S = energy slope which, for the case of quasi-uniform flow, equals the bed slope S b and water surface slope S w , and z is the vertical coordinate.Equation 3 is considered valid at an elevation away enough from the bed roughness crest z c .Note that, if C d and d p are vertically uniform, 〈u〉 IL is constant over depth.
Another important flow feature in vegetation canopies is the turbulent kinetic energy, defined as , where u′, v′, and w′ are the temporal velocity fluctuations in streamwise, spanwise and vertical directions respectively.Tanino and Nepf (2008) proposed a model to predict the turbulent kinetic energy in the IL 〈k t 〉 IL assuming that the main source of turbulence production in this layer is associated with the vortex shedding from individual cylinders: where γ ≈ 1.1 is an empirical constant for sparse arrays (Tanino & Nepf, 2008).
In the NBL, the no-slip condition and bed roughness have a significant influence on the flow structure leading to significant vertical gradients of streamwise velocity.Higher ratios of streamwise velocity and bulk flow velocity (〈u〉/U) were observed close to the bed roughness crest in vegetated regions compared to non-vegetated regions for the same discharge and bed slope (Penna et al., 2020).In addition, the vertical distribution of the doubleaveraged streamwise velocity in the NBL resembles a logarithmic profile for both smooth bed (Conde-Frias Water Resources Research 10. 1029/2023WR035879 et al., 2023) ) and rough beds (Penna et al., 2020;Ricardo et al., 2016).Conde-Frias et al. (2023) defined the thickness of the NBL δ as the elevation where the streamwise velocity reaches the vertically uniform velocity, that is, 〈u〉 z=δ = 〈u〉 IL .For vegetated smooth bed flows, δ is negatively correlated with stem array density (Conde-Frias et al., 2023) and is therefore negatively correlated with 〈k t 〉 IL .
The NBL turbulence structure in vegetated rough bed flows is controlled by the bed roughness and vegetation (Penna et al., 2020).Peak values of the spatially-averaged turbulent kinetic energy and turbulent stresses have been found close to the bed, and these quantities decrease up to the IL where they become constant (Penna et al., 2020;Ricardo et al., 2016).In addition, the streamwise and spanwise fluctuations of turbulent kinetic energy in the NBL are small compared to the spatial fluctuation of turbulent kinetic energy in the IL, suggesting different types of vortices governing velocity fluctuation in each layer (Penna et al., 2020).Furthermore, spectral analysis showed that macroscale turbulence (characterized by low-wave numbers) close to the bed does not change between vegetated and unvegetated regions and depends only on the grain median diameter (Penna et al., 2020).Yang et al. (2016) proposed a model assuming that the total turbulent kinetic energy 〈k t 〉 NBL in the NBL corresponds to the sum of bed-generated turbulence 〈k t,b 〉 and vegetation-generated turbulence 〈k t,v 〉 as They recommended estimating 〈k t,b 〉 using empirical relationships relating turbulent kinetic energy to bulk flow velocity for bare bed conditions.For the estimation of 〈k t,v 〉, Yang et al. (2016) assumed that 〈k t,v 〉 ≈ 〈k t 〉 IL .

Linkages Between Bed Shear Stress and Flow Characteristics
Bed shear stress in vegetated areas can be analytically assessed considering a force balance, that is, the so-called depth-slope method (e.g., Aberle & Järvelä, 2013).For the case of steady quasi-uniform flow for plane beds (in the absence of bedforms), bed shear stress τ b can be expressed as where ϕ b = solid volume fraction occupied by the bed, z t = elevation of roughness troughs, and 〈h〉 = 〈z ws z c 〉 (See Figure 1).The first term on the right-hand side of the equation represents the total shear stress (τ t ), and the second term represents the vegetative drag force per unit bed area.Although this method allows the estimation of bed shear stress when the total shear stress and drag forces are known (Armanini & Cavedon, 2019;Jordanova & James, 2003;Wang et al., 2023), it does not provide insight into how bed shear stress relates with the changes in the flow structure.
Following DAM, bed shear stress can be expressed in terms of the near-bed flow field assuming that the momentum balance close to the bed roughness crest is governed by the turbulent and form-induced stresses (Nikora et al., 2007) according to; This method relies on an accurate characterization of the flow field in the NBL.Therefore, it requires measurements with high spatial resolution inside the canopy (Xavier & Wilson, 2007).
Moreover, there is a strong correlation between bed shear stress and squared bulk velocity.A common approach developed for unobstructed open-channel flows and adopted to vegetated flow (e.g., Baptist et al., 2007;Le Bouteiller & Venditti, 2015) is to scale bed shear stress with the squared bulk flow velocity through empirical coefficients such as the Darcy-Weisbach bed friction factor λ b as Water Resources Research

10.1029/2023WR035879
For open-channel flows, λ b depends on the Reynolds number and the ratio between bed roughness height and water depth (e.g., Yen, 2002).Although Schoneboom et al. (2010) showed that the presence of vegetation affects λ b , there are no expressions available to date to directly estimate λ b in vegetated flows.
Alternatively, bed shear stress can be expressed in terms of near-bed turbulent kinetic energy (Kim et al., 2000) as: where the scaling constant C k t = 0.19 was derived for tidal flows but works also for open-channel flows (Biron et al., 2004;Rowiński et al., 2005).However, considering that turbulent kinetic energy in vegetated flows is not only generated by the bed friction, the presence of vegetation may have a significant impact on C k t .
Based on the characterization of the flow structure in the NBL in vegetated smooth beds, Conde-Frias et al. ( 2023) recommended determining bed shear stress from where μ = dynamic viscosity, and C δ = scaling coefficient (1.98 for ϕ v > 0.04).This model is valid for smooth beds and it is based on the correlation between viscous stresses and the ratio between 〈u〉 IL and δ.Although this model was derived for smooth beds, a proper characterization of the combined impact of vegetation and bed roughness on δ can provide insight into the description of bed shear stress in vegetated rough bed flows.

Experimental Setup
Experiments were carried out in a 32-m long, 0.4-m deep, and 0.6-m wide tilting flume at the Leichtweiß-Institute for Hydraulic Engineering, Technische Universität Braunschweig, Germany (Figure 2).The flow discharge Q, supplied by the laboratory water circuit, was measured in the supply pipe through an inductive flow meter (Krohne Aquaflux 090K, 0.3% accuracy).The water depth was regulated by a tailgate, located 24.4 m downstream of the flume inlet.The water surface elevation z ws was measured through a system of 12 piezometers installed at the flume bottom at distances of 12, 13, 14, 15, 15.15, 15.35, 15.55, 15.95, 16, 16.15, 17, and 18 m from the flume inlet, through communicating pipes.A movable ultrasonic sensor (Microsonic mic + 130/IU/TC, 0.15% accuracy) measured the water elevation in each communicating pipe for 30 s.The measurements were repeated every 6 min until the measurements converged to a constant value (±1 mm).This high spatial and temporal resolution allowed for a proper characterization of the water surface considering that, especially in emergent conditions, vegetated flow is characterized by strong undulating water surface.The water surface slope S w was estimated through linear regression of the measured water surface elevations using MATLAB's fitlm.m function.
Rough bed conditions were exemplarily simulated by covering the flume bottom with an artificial grass mat, which is a common approach to represent understory grasses commonly found in floodplains (Box et al., 2021).The 16 m long mat was installed at a distance of 7 m from the flume inlet and consisted of 0.026 m high blades with an array density of approximately 492,800 stems/m 2 and a solid volume fraction ϕ b occupied by the grass of 0.059 (Branß & Aberle, 2022).This grass mat, which served as bed roughness in the present study, represents a dense grass canopy according to Nepf (2012).Only one bed roughness was used since the focus of the study was on the evaluation of the impact of the array density on bed shear stress.
Rigid emergent vegetation was mimicked by 0.3 m high plastic rigid cylinders with a diameter of d v = 0.01 m.
The cylinder array covered a distance of 8 m starting at a distance of 12 m from the flume inlet.Following Lindner (1982), this canopy length suggests that at the middle of the used canopy length, the hydraulic conditions are equivalent to a canopy of infinite length.The cylinders were arranged in a staggered pattern to have closer hydrodynamic behavior to random arrays found in nature.Four canopy configurations were tested by varying the cylinder distance ∆s (0.40, 0.30, 0.20, and 0.15 m) resulting in stem areal densities m of 6. 25, 11.11, 25, and 44.44 stems/m 2 (Figure 3).All configurations can be classified as sparse vegetation since d v / Δs < 0.56.
A total of five experimental series were carried out using the aforementioned four canopy configurations and a reference series without vegetation.In each series, five flow discharges Q were imposed ranging between 0.03 and 0.09 m 3 /s using a constant spatially averaged flow depth 〈h〉 of approximately 0.224 m, defined as the distance from the mean water surface z ws to the top of the grass mat height z c (cf. Figure 1).Thus, the distance from the mean water surface z ws to the flume bottom z t was 0.25 m.These hydraulic conditions were chosen because they are comparable to measurements carried out by Schoneboom et al. (2011) in the same flume and with the same emergent rigid cylinders with an array density of 25 stems/m 2 , discharges between Q = 0.03-0.105m 3 /s, and a water depth of z ws z t = 0.25 m.Quasi-uniform flow conditions were achieved by adjusting the bed slope and tailgate iteratively for each experiment ensuring that the bed slope matched the water surface slope.Stem Reynolds numbers were in the range for which C d , is constant (Schoneboom et al., 2011).Finally, the Froude number Fr = U/ ̅̅̅̅̅̅̅̅̅ g〈h〉 √ was always well below 1 so that the flow was always subcritical.Table 1 summarizes the hydraulic boundary conditions.For each array density, velocity measurements in the near-bed region were carried out in 4 bed-parallel planes that were located between 2.5 and 12.5 mm above z c , leading to elevations relative to the water depth of (z z c )/ 〈h〉 = 0.01-0.06.These were chosen according to ranges reported in previous research dealing with the NBL in flows with rigid emergent vegetation and rough beds, (z z c )/〈h〉 = 0.005-0.2(Penna et al., 2020;Yang et al., 2016;Yang & Nepf, 2019).In each plane, 35 velocity time series were recorded at points defined by a square grid with streamwise (Δx) and spanwise (Δy) resolutions of 150 mm, and 60 mm (Figure 3).The velocity measurements were carried out in coincident mode up to 20,000 samples or 90 s.Data records with less than 10,000 samples were discarded, leading to data records with (on average) 17,000 samples and sampling frequencies of 270 Hz.
The spatial resolution of the measurements was chosen by comparing the spatially averaged streamwise velocity and turbulent flow quantities 〈k t 〉, ⟨u′w′⟩ , and 〈 ũ w〉 to values obtained with a denser reference staggered grid (Δx = 100 mm, Δy = 30 mm) of 62 vertical profiles following Xavier and Wilson (2007) using a discharge of 0.06 m 3 /s and 6.25 stems/m 2 .The absolute relative errors between the sparse and dense grid were of 0.00%, 4.0%, 3.4%, and 8.1% for 〈u〉 , 〈k t 〉, ⟨u′w′⟩ , and 〈 ũ w〉 , so it was decided to use the sparse grid for the measurements (further details can be found in Text S1 in Supporting Information S1).Finally, temporal convergence was validated by bootstrapping analysis (Buffin-Bélanger & Roy, 2005), showing that 〈u〉 , 〈k t 〉, and ⟨u′w′⟩ , converged (±5%) after 10,000 samples.
Temporal post-processing of the velocity series was carried out using the LDA Toolbox of Aleixo et al. (2016).Spatially averaged quantities were calculated by averaging all temporally post-processed samples for each plane parallel to the bed.

Bed Shear Stress Measurements
Bed shear stress was estimated through a shear plate by measuring the streamwise hydrodynamic forces acting on the bed, F sp , through a force sensor connected to a freely moving plate.Although this method has been successfully applied in other hydraulic fields (Park et al., 2019;Pujara & Liu, 2014), its application in vegetated flows is scarce (Ishikawa et al., 2003;Tinoco & Cowen, 2013).The used shear plate is an adaptation of a previous design of Niewerth et al. (2021) and consists of a movable wooden plate connected to the top of a drag force sensor (DFS), as shown in Figure 4.The DFS was developed by Schoneboom et al. (2008) and has previously been used to measure the drag force exerted by vegetation elements in various studies (Jalonen et al., 2013;Siniscalchi et al., 2012;Sukhodolov et al., 2022;Västilä et al., 2013;Vettori et al., 2021).In brief, it estimates forces by measuring the deformation of a 3 mm thick bending stainless-steel beam through eight strain gauges configured as two Wheatstone bridge configurations.More details on the DFS can be found in the above references.
The 0.46 m long and 0.36 m wide shear plate was located at a distance of 15.6 m from the flume inlet in the middle of the flume.A 5 mm wide gap between the shear plate and the channel bottom ensured free movement of the shear plate.Cylinders were mounted onto the shear plate, allowing for the measurement of the combined effects of the drag forces acting on the cylinders attached to the shear plate and the bed shear stress acting on the surface of the plate.For the case of 25 stems/m 2 , only one cylinder was located on the shear plate because the remaining four cylinders were too close to the edges of the shear plate affecting its stability (cf. Figure 3).Force measurements were carried out with a sampling frequency of 1,613 Hz for 180 s and repeated three times in 45 min.
Bed shear stress (τ b,sp ) was estimated by subtracting the drag forces acting on the cylinders attached to the shear plate from the drag force measured by the shear plate and dividing this difference by the shear plate area, that is, where n sp = number of cylinders mounted on the shear plate and A sp = area of the shear plate.Analogously, the total shear stress (τ t,sp ) can be estimated by where the second term on the right side of the equation corrects the difference between the array density in the flume and the shear plate since the drag forces acting on the cylinders on the shear plate are already considered by F sp .
For the determination of the drag forces, the Lindner (1982) method was used to estimate C d .This approach was chosen since Schoneboom et al. (2011), using the same cylinders and flume, found that this approach yields good results compared to direct drag force measurements, We therefore used C d = 0.94, 0.98, 1.07, and 1.17 for the array densities of 6. 25, 11.11, 25, 44.44 stems/m 2 , respectively, to estimate 〈F d 〉 according to Equation 1.
To assess the performance of the shear plate in vegetated flows, additional estimations of bed shear stress (τ b,hs ) and total shear stress (τ t,hs ) were obtained using the depth-slope method according to Equation 6, which relies only on a momentum balance in contrast to the other aforementioned methods (cf.Section 2.3).This assessment was done using uncertainty analysis (UA) following the "Guide of Expression of Uncertainty in Measurement" (GUM) (JCGM, 2008), which has been recently recommended by the hydraulic community (Mrokowska et al., 2013;Muste et al., 2017).The combined uncertainty of a measurand Y = f(X i , X i+1 , …, X n ) can be estimated as where Y is the measurand, X i are the dependent variables, and s denotes the variable uncertainty.The combined uncertainty of bed shear stress using the shear plate (s(τ t,sp )) is estimated as where the force measuring uncertainty was assessed by applying a set of known loads through a pulley system at static water conditions (Q = 0 m 3 /s, S b = 0, 〈h〉 = 0.224 m) leading to an uncertainty of s(F sp ) = 0.03 N (see Text S2 in Supporting Information S1).The force measurement accuracy could also be affected by flowing water due to the presence of flow gradients at the edges of the shear plate (Pujara & Liu, 2014).Therefore, velocity measurements were carried out to determine the presence of streamwise flow gradients induced by the gaps at the edges of shear plates.The measurements showed no significant spatial fluctuations of u(z) and k t (z) above the gap compared to the spatial fluctuations in the streamwise direction along the centerline of the shear plate (see Text S3 in Supporting Information S1).Water Resources Research 10.1029/2023WR035879 Analogously, the combined uncertainty of bed shear stress using the depth-slope method (s(τ b,hs )) is The uncertainty related to the water slope was estimated from the 95% confidence interval of the slope obtained from linear regression.Therefore the uncertainty of τ b,hs is also dependent on the flow spatial heterogeneity through all its sources of uncertainty.This means that higher uncertainties are expected using the depth-slope method at the flow conditions where the spatial flow heterogeneity is higher.

Validation of Bed Shear Stress Estimations
Quasi-uniform flow was achieved for all runs (Figure 5a), with a mean absolute percentual error (MAPE) between water surface slopes and bed slopes of 7.5%.We acknowledge that the UA confirmed that, besides our effort to provide a higher spatial resolution of z ws , the water surface slope is prone to errors in emergent flow conditions due to the undulating water surface (relative uncertainties between 12% and 36%).On the other hand, for all runs the desired mean water depth of 0.224 m was achieved with low uncertainties (0.2%-4.5%).
We observed a general disagreement between the shear plate (τ b,sp ) and the depth-slope (τ b,hs ) method for the highest array density (Figure 5b).As shown in Table 2, MAPE ranges between 10.9% and 167.4%.The UA revealed that bed shear stress estimation using the shear plate is more precise than using the depth-slope method (Table 2), with less uncertainty both in absolute and relative terms.The uncertainties in the depth-slope method can bed mainly attributed to the uncertainties associated with the water surface slope estimation, whereas the uncertainty of the shear plate method is mainly attributed to the uncertainty associated with the drag force estimation.Due to higher precision of the shear plate, bed shear stress and total shear stress estimations using the shear plate are used, and referred to as τ b and τ t respectively.

Impact of Vegetation on Bed Shear Stress
For constant water depth and discharges ranging between 0.03 and 0.09 m 3 /s, increasing array density led to higher required bed slopes and consequently water surfaces slope to achieve quasi-uniform flow (Figure 6a).We acknowledge that this does not reflect natural conditions (constant bed slope)  for which the increase of array density will lead to an increase in water depth.Nonetheless, the chosen measurement strategy allows to study the impact of vegetation for a range of bulk flow velocities.Consequently, for comparable bulk flow velocities, the total shear stress increased with array density, due to higher resistance of vegetation per unit bed area (〈F d 〉 m), which in turn led to higher required water surface slopes (Figure 6b).
In addition, Figures 6a and 6b show that the acquired data are consistent with the data from Schoneboom et al. (2011) who used similar hydraulic conditions (cf.Section 3.1) and the same cylinder with an array density of 25 stems/m 2 (with in-line [Sch-L] and staggered [Sch-S] arrangement).However, compared to the present experiments, the experiments of Schoneboom et al. (2011) were carried out with a smoother bed composed of 3 mm high rubber pyramids and a water depth of 0.25 m.This depth is 12% higher than in the present experiments when the roughness crest z c is taken as the reference level, but is the same when the roughness troughs z t define the reference level.The experiments of et al. ( 2011) with cylinder arrays in a staggered pattern required a higher water surface slope to achieve quasi-uniform flow compared to the in-line pattern (Figure 6a).This was due to higher drag exerted by cylinders in a staggered pattern which, in combination with the higher slope, lead to higher total shear stress (Figure 6b).
As shown in Figure 7a, which presents the ratio between bed shear stress and total shear stress τ b /τ t as a function of bulk velocity U, an increase in array density lead to a decrease in τ b /τ t .This in turn means the contribution of bed shear stress to total shear stress is not neglectable for the investigated sparse arrays of vegetation, which is consistent with the results of Schoneboom et al. (2011).Note that, although the experiments of Schoneboom et al. (2011) were carried out with 25 stems/m 2 , the lower τ b /τ t values compared to the series m25 can be attributed to the fact the bed roughness was smoother and the water depth above the roughness crest was higher.Figure 7b shows that the increase in array density led to higher bed shear stress for a constant water depth and comparable bulk flow velocities.Similar to unobstructed open-channel flows, there is a strong correlation between bed shear stress and the square of bulk flow velocity, with Pearson's correlation coefficient (r) of 0.99 for all series, as well as for the Schoneboom et al. ( 2011) series (r = 0.97).Note that correlation coefficients are labeled as weak (|r| ≤ 0.5), moderate (0.5 < |r| < 0.8), and strong (|r| ≥ 0.8) following Devore (2012).It needs to be mentioned again that the experiments were carried out at different bed slopes that the bed shear stress results must be interpreted in light of different total shear stress.However, there is a high variation in the average increment of bed shear stress between consecutive series with increments of bed shear stress of 60%, 1%, 21%, and 86% between the series with 0, 6. 25, 11.11, 25, and 44.44 stems/m 2 .Note that a repetition of the experiments after a detailed inspection of the shear plate led to the same results, so that these variations are not attributed to measurement errors.
We therefore infer that the data reveal an increase in the ratio of bed shear stress and bulk flow velocity with increasing array density which can be further investigated based on the Darcy-Weisbach bed friction factor (λ b ).Figures 8a and 8b show λ b for the different series (cf.Equation 8) dependent on discharge and slope.We found that λ b increases with array density (r = 0.95), and, consistent with Schoneboom et al. (2010), that the bed friction factor is higher for vegetated flows compared to non-vegetated flows.
There is also a strong correlation of λ b with flow discharge (r = 0.81-0.96)(and hence bulk velocity) which is consistent with measurements of Box et al. (2021) for the understory grasses.Box et al. (2021) attributed this increase in bed friction factor to a higher momentum penetration in the grass with higher velocities, which is equivalent to an increase in effective bottom roughness.
However, there is also a strong correlation of λ b with bed slope (r = 0.74-0.94)for each series, except for the series m25 (r = 0.09) suggesting that increases of λ b can also be related to changes in bed slope.Notwithstanding, because in our experiments change in bed slope were required to keep a constant submergence ratio between water depth and roughness height for each array density, we cannot state clearly if there is a strong dependency of bed slope over λ b independent of array density.
However, when all data are correlated with hydraulic conditions, lower correlations are found (r = 0.02 and 0.55 for discharge and bed slope, respectively), showing that the changes in λ b are governed by the array density.In this context, Figure 8b shows that for each bed slope there is a clear increase in λ b with array density.
Thus, the results suggest competing impacts of vegetation on bed shear stress regarding the investigated constant water depth scenario.On one hand, the ratio between bed shear stress and bulk flow velocity increases with array density being the reason for the increase of the bed friction factor.On the other hand, in natural conditions with a fixed slope, the presence of vegetation with increasing array densities leads to an increase in water depth and therefore a reduction of bulk flow velocities.Hence it may be expected that both bed shear stress and bed friction factors decrease.Therefore, although the presence of vegetation may lead to a net reduction of bed shear stress, the effect of array density on the bed friction factor hampers the estimation of the total reduction of bed shear stress.
This conclusion is consistent with the bed shear stress model presented by Lu et al. (2021) for rigid emergent vegetation, in which the impact of array density on bed shear stress was analytically considered.This model is based on the phenomenological theory turbulence, which states that the ratio between bed shear stress and bulk flow velocity is governed by the size ratio between macroscale eddies to sediment size eddies due to the energy cascade process.Lu et al. (2021) argue that the presence of vegetation constrains the size of the macroscale eddies so that the ratio between macroscale eddies and sediment size eddies decreases with array density, leading to a higher energy transfer between the bulk flow velocity to the near-bed velocity governing bed shear stress.Therefore, although vegetation leads to a decrease in bulk flow velocity, it leads to an increase in the ratio of bed shear stress to bulk flow velocity which is proportional to array density.This is also confirmed by our findings.

Streamwise Velocity
Vertical profiles of near-bed double-averaged streamwise velocity normalized by bulk flow velocity 〈u〉/ U are presented for all array densities in Figure 9.The figure shows that, for each series, 〈u〉/U resembles a logarithmic profile which is consistent with previous observations in experiments with vegetated rough beds (Penna et al., 2020;Ricardo et al., 2016).The logarithmic profile is described as where C 1 , C 2 , and C 3 are fitting coefficients.All series present high coefficients of determination (R 2 ) (Table 3), except the m06 series (6.25 stems/m 2 ) which is characterized by a slightly lower R 2 -value.In addition, 〈u〉/U increases with array density in the NBL, consistent with results presented by Conde-Frias et al. (2023) for the case of vegetated smooth beds.
Based on the logarithmic fits, the NBL thickness (δ) can be estimated for each array density.Following Conde-Frias et al. (2023), δ corresponds to the first elevation above z c where the streamwise velocity reaches the magnitude of bulk flow velocity, that is, δ = z z c when 〈u〉 z z c /U = 1.Here, the relative NBL thickness (δ/〈h〉) can be estimated using Equation 16 at the boundary condition 〈u〉 (z z c ) / 〈h〉 / U = 1 from Table 3 reports the estimated values of δ/〈h〉 and these results show that for sparse arrays of vegetation, the thickness of the NBL can represent an important portion of the water depth, which has not been previously characterized for vegetated rough bed flows.As shown in Figure 10a, δ/h is inversely proportional to array density (r = 0.84), which is consistent with the results of Conde-Frias et al. ( 2023) for smooth beds.In addition, there is a strong correlation between δ/h and the estimations of the vegetation-induced turbulence 〈k t,v 〉 (r = 0.86) (Figure 10b).Therefore, as stated by Conde-Frias et al. (2023), the thickness of the NBL seems to be controlled by the vertical momentum transfer induced by the turbulent kinetic energy production in the IL.
It is observed that the NBL thickness δ decreases with array density, and considering that our results shows an increase of bed shear stress with array density for constant U and 〈h〉, it can be assumed that bed shear stress scales with 1/δ.This is consistent not only with empirical findings of Conde-Frias et al. (2023) for smooth beds but also with the Boundary Layer Theory (e.g., Schlichting & Gersten, 2017) stating that bed shear stress scales not only with  the square of the velocity in the undisturbed region, but also decreases with the thickness of the boundary layer.

Turbulent Kinetic Energy
The vertical distribution of normalized values of turbulent kinetic energy 〈k t 〉/U 2 is shown in Figure 11.The data represent mean values for each array density for all discharge conditions, as there is no strong and consistent correlation between these quantities and flow discharge (Table 4).We note that in regard to the vertical gradient, there is only a moderate to strong correlation observable for the m25 and m44, whereas, for the series m06 and m11, there is a weak correlation with elevation.Therefore, the near-bed value of turbulent kinetic energy 〈k t 〉 NBL / U 2 for each series is defined as the average of 〈k t 〉/U 2 for all discharges at the closest elevation to the bed, that is, (z z c )/〈h〉 = 1.3%.Furthermore, there is a weak impact of array density on turbulent kinetic energy (r = 0.47), and its impact is in the order of magnitude of the variation due to flow discharge (error bars) and elevation.
In addition, measured values of 〈k t 〉 NBL / U 2 were compared against predicted values the near-bed turbulent kinetic energy model of Yang et al. (2016) (Equation 5) (Figure 12).The bed roughness-induced turbulence 〈k t,b 〉 was estimated through the scaling relation between turbulent kinetic energy and bed shear stress for open channel flows (Equation 9) using the bed shear stress values for the bare bed conditions.The vegetation-induced turbulence 〈k t,v 〉 was estimated as the turbulent kinetic energy in the IL (Equation 4). Figure 12 shows that the 〈k t 〉 NBL predictions are in the same order of magnitude as the measured values.The Yang et al. (2016) model results in larger values than the measurements (MAPE = 52%), and this overestimation can be attributed to the prediction of 〈k t,b 〉.Furthermore, both the measurements and the predicted values show that the variation of 〈k t 〉 NBL with array density is small compared to the magnitude of the bed roughness-induced turbulence, with mean ratios 〈k t,v 〉/〈k t,b 〉 of 3.1, 4.7, 8.5, 13.2% for the array densities of 6.25, 11.11, 25, 44 stems/m 2 respectively.Thus, consistent with the measurements of Penna et al. (2020), and the predictions of Yang and Nepf (2019), it can be concluded that near-bed turbulent kinetic energy is dominated by the bed roughness-induced turbulence for sparse arrays of rigid emergent vegetation.

Turbulent and Form-Induced Stresses
Vertical distributions of the normalized streamwise turbulent stresses ⟨u′w′⟩/ U 2 and form-induced stresses 〈 ũ w〉/ U 2 are shown in Figure 13.Similarly to turbulent kinetic energy, there is a weak correlation between turbulent stresses as well as form-induced stresses and discharge (Table 5).In regard to the vertical distribution in the near-bed region, the data do not show a consistent correlation as observed for the turbulent kinetic energy.Therefore, in analogy to the turbulent kinetic energy, near-bed values were defined as the average for all flow discharges at the closest elevation to the bed.Regarding the impact of array density, there is a weak correlation between array density and turbulent stresses (r = 0.11), suggesting that turbulent stresses are controlled by bed roughness.However, there is a moderate impact of array density on form-induced stresses (r = 0.56).In addition, forminduced stresses are significantly smaller than turbulent stresses (on average 7% of turbulent stresses).Although previous measurements in vegetated flows show comparable values of turbulent stress and form-induced stress close to the bed (Ricardo et al., 2016) (m = 400-1,600 stems/m 2 , mixture of sand and gravel), it seems that for sparse arrays of vegetation, the array density is too small to generate form-induced stresses comparable to turbulent stresses.

Linkages Between Near-Bed Flow and Bed Shear Stress
To analyze which flow characteristics relate better with bed shear stress independent of array density, bed shear stress was plotted against different flow variables 14).First, as a reference case, bed shear stress was correlated with the square of bulk flow velocity (Figure 14a).For this case, a strong correlation is found (r = 0.84), but higher bed shear stress values are observed for higher array densities, especially for the case of 44.44 stems/m 2 .As mentioned in Section 4.1, these higher bed shear stress values are not resulting from measurement errors.
Following Kim et al. (2000) (Equation 9), bed shear stress is compared to near-bed values of turbulent kinetic energy in Figure 14b.The correlation between near-bed turbulent kinetic energy and bed shear stress is similar to the square of bulk flow velocity (r = 0.86).This is attributed to the fact that there is no significant impact of vegetation on the near-bed turbulent kinetic energy.On the other hand, when bed shear stress is compared to the vegetation-induced component of turbulent kinetic energy (Figure 14c), there is a better agreement with r = 0.91.This suggests that the turbulent production at the IL impacts bed shear stress due to higher vertical momentum exchanges as pointed out by Conde-Frias et al. (2023).
The sum of near-bed values of turbulent and form-induced stress is compared to bed shear stress following Equation 7 in Figure 14d.Although this sum gives values of the same order of magnitude, as turbulent stress is not affected by array density, the correlation is similar to the previous quantities (r = 0.89).
Moreover, when bed shear stress is compared to the ratio between the squared bulk flow velocity and the relative thickness of the NBL, a better agreement is observed (r = 0.93) leading to a more unified trend between all series (Figure 14e).This is attributed to the fact that δ/h serves as a scaling parameter between bulk flow velocities and near-bed velocity, as was previously observed for the case of vegetated smooth beds (Conde-Frias et al., 2023).Assuming that bed shear stress is governed by the streamwise velocity acting just above the bed roughness, the  relative thickness of the NBL, therefore, serves as a scaling parameter to relate bed shear stress with the square of bulk flow velocity.

Conclusions
In this paper, we explored the impact of array density on bed shear stress and near-bed flow structure in rough beds covered with sparse arrays of rigid emergent vegetation.Through the combination of novel direct measurements of bed shear stress and detailed characterization of the near-bed flow layer, we quantified how the relationships between bed shear stress and flow field are affected by the presence of vegetation in rough beds, and how bed shear stress relates to near-bed flow quantities, highlighting how a small concentration of elements can significantly modify the bed shear stress and near-bed flow field.
Regarding the impact of vegetation on bed shear stress and bulk flow characteristics, an increase in the array density enhances the ratio between bed shear stress and bulk flow velocity for constant water depths, leading to higher values of Darcy-Weisbach bed friction factor.This behavior suggests that the presence of vegetation increases the energy transfer toward the bed.However, further research is required to analyze the impact of array density on bed shear stress at constant slope conditions, that is, increase of water depth with array density.
With respect to the changes in the near-bed flow structure, higher ratios of near-bed streamwise velocity to bulk flow velocity are observed for higher array densities.The results show that the thickness of the NBL decreases with array density, as occurs for smooth bed flows.In addition, it was observed that for sparse arrays near-bed turbulence is governed by the bed roughness with minor contributions of vegetation-induced turbulence in the NBL.However, changes in the NBL flow structures seem to be significantly influenced by the turbulence production at the IL.
When comparing bed shear stress to different flow descriptors, it was found that the NBL thickness is a proper characteristic length to scale bed shear stress with the bulk flow velocity.Therefore, we conclude that bed shear stress in vegetated rough beds is governed by near-bed streamwise velocities rather than near-bed turbulent descriptors.However, turbulence production at the IL, which is affected by array density, seems to have an impact on the near-bed streamwise velocity and hence bed shear stress.Further research will explore in detail the linkages between turbulence production over the full water depth in a wider range of hydraulic conditions to enhance the understanding of the impact of vegetation on the NBL flow structure and bed shear stress.

List of Symbols
A sp area of the shear plate

Figure 1 .
Figure 1.Schematic representation of the flow field in emergent rigid vegetation flows.Dark blue and red lines represent the vertical distribution of the double-averaged streamwise velocity 〈u〉 and turbulent kinetic energy 〈k t 〉 in contrast to open-channel flow conditions (light lines).

Figure 2 .
Figure 2. Side view of the main components and measurement system at the tilting flume (not to scale).

Figure 3 .
Figure 3. Cylinder arrangement for all array densities and location of vertical profiles (crosses).The dotted line represents the location of the shear plate (cf.Section 3.3).For the array m25, four stems were moved out of the shear plate to increase the stability of the shear plate (empty circles).

Figure 4 .
Figure 4. Working principle and components of the shear plate.
The uncertainty in drag force s(〈F d 〉) is related to the uncertainty associated with C d and 〈h〉.The uncertainty of C d was estimated based on the standard deviation of the spatial fluctuation of the 〈F d 〉 measured by Schoneboom et al. (2011), leading to s C d ) = 0.08C d .The uncertainty of 〈h〉 was estimated according to the standard deviation of the spatial fluctuations of h.The uncertainty of 〈F d 〉, and consequently τ b,sp , depends on the spatial heterogeneity of the flow depth governing the frontal area exposed to the flow which in turn is required for the estimation of 〈F d 〉.

Figure 5 .
Figure 5. Agreement between bed slope and water surface slope (a).Agreement and uncertainty of bed shear stress estimations through the depth-slope (τ b,hs ) and the shear plate method (τ b,sp ) (b).Dotted lines denote perfect agreement.Note that bed shear stress estimations for m44 at Q = 0.09 m 3 /s were discarded because the shear plate reached the maximal displacement.

Figure 6 .
Figure 6.Required water surface slope for quasi-uniform flows for different discharges and array density (a).Impact of array density over total shear stress (b).Same water depth for all experiments.Absence of an error bar means that the uncertainty is smaller than symbol size.

Figure 7 .
Figure 7. Impact of stem areal density on the ratio bed shear stress to total shear stress (a) and on bed shear stress (b).Error bars in (b) represent the uncertainty of τ b .Absence of an error bar means that the uncertainty is smaller than the symbol size.

Figure 8 .
Figure 8. Impact of array density on Darcy-Weisbach bed friction coefficient λ b in terms of flow discharge (a), and bed slope (b).Error bars represent the combined uncertainty of λ b due to uncertainty propagation according to Equation 13.

Figure 9 .
Figure 9. Logarithmic fits of vertical profiles of streamwise velocity for all array densities.Data points represent the average for all flow discharges.Shaded areas represent the uncertainties based on the root mean square error (RMSE) of each fit.

Figure 10 .
Figure 10.Impact of array density (a) and vegetation-induced turbulent kinetic energy (b) on near-bed layer relative thickness.Error bars represent the uncertainty of δ/〈h〉 based on the root mean square error (RMSE) of each fit.

Figure 11 .
Figure 11.Vertical profiles of normalized turbulent kinetic energy for all series.Error bars represent the standard deviation due to different flow discharges.

Figure 12 .
Figure 12.Comparison of measured versus predicted spatially averaged near-bed turbulent kinetic energy (Equation 5).Error bars of measured 〈k t 〉 NBL represent the standard deviation due to different flow discharges, whereas for estimated 〈k t 〉 NBL represent the propagated uncertainty in the model.

Figure 13 .
Figure 13.Vertical distribution of streamwise components of turbulent stresses (a), and form-induced stresses (b).
an isolated cylindrical object C d mean bulk drag coefficient of the array C k t scaling constant for τ b prediction in terms of k t (Equation9)C δ scaling constant for τ b prediction in terms of δ (Equation10)C 1 , C 2 ,C 3 fitting coefficients of the logarithmic profile (acting on the bed f b sum of viscous and pressure drag per unit fluid volume due to bed roughness f v sum of viscous and pressure drag per unit fluid volume due to vegetation g i i-to component of the gravity acceleration h flow depth

Figure 14 .
Figure 14.Correlations between bed shear stress and flow characteristics.(a) Bulk flow velocity, (b) near-bed turbulent kinetic energy, (c) vegetation-induced turbulent kinetic energy, (d) near-bed Reynolds and form-induced stresses, (e) ratio between bulk flow to near-bed layer thickness.

Table 2
Quality Indicators of Bed Shear Stress Estimations Note.MAPE, mean absolute percentage error between τ b,sp and τ b,hs ; s, measurand uncertainty.

Table 3
Fitting Coefficients of Logarithmic Fit for All Series

Table 4
Correlation Coefficients Between Turbulent Kinetic Energy With Flow Discharge and Relative Elevation