3.2.1.1. Effects of varying BC, d₅₀/h, VR and the mussels’ arrangement on the bed surface

The vertical distributions of the double-averaged streamwise velocity, $\langle \overline{u}\rangle$ , are qualitatively similar in the mussel-bed simulations performed with BC ≤ 0.280 and with various values of VR and $ d_{50}/h $ (figure 7). These profiles are characterised by a rapid increase of $\langle \overline{u}\rangle$ with the distance from the location where $\langle \overline{u}\rangle$ = 0. The rate of increase of $\langle \overline{u}\rangle$ inside the inner layer decreases with increasing surface mussel coverage density (figure 7 a) because the mussel-induced drag force per unit area is proportional to BC. By contrast, once the skimming flow regime becomes dominant, $\langle \overline{u}\rangle$ ≈ 0 over a significant fraction of the inner layer (e.g. until $ z/h $ ≈ 0.3 for the simulation conducted with BC = 0.828 shown in figure 7 b). Above the level where $\langle \overline{u}\rangle$ ≈ 0, the rate of increase of $\langle \overline{u}\rangle$ inside the inner layer in the BC = 0.828 simulation is comparable to that observed in the simulation performed in a channel with no mussels and $ d_{50}/h $ = 0.13 (figure 7 a). This suggests that the skimming flow regime becomes dominant when the packing of the mussels is so high that mussels can be considered to behave similarly to standard rough surfaces with uniformly distributed roughness. Given that the channel discharge was the same in all simulations, a decrease of $\langle \overline{u}\rangle$ inside the inner layer is accompanied by an increase of $\langle \overline{u}\rangle$ inside the top part of the inertial layer. For cases with sparsely distributed mussels, the level where the flow inside the inertial layer starts accelerating due to increasing BC is situated at z ≈ 2h.

The presence of distributed bed roughness increases the drag acting on the bottom fluid inside the inner layer, which, in turn, decreases the rate of increase of $\langle \overline{u}\rangle$ inside the inner layer (figure 7). As BC increases and the transition to the skimming flow regime starts, this effect weakens as most of the gravels are situated inside a region containing slowly moving fluid. For example, the effect of replacing the smooth channel bed by a gravel surface with $ d_{50}/h $ = 0.13 becomes negligible in the BC = 0.280 simulations (figure 7 c).

The effect of the active filtering on the vertical profile of $\langle \overline{u}\rangle$ is relatively minor as exemplified in figure 7(d), which compares the BC = 0.280, $ d_{50}/h $ = 0.13 simulations conduced with VR = 0.5 and VR = 1.22. In the presence of strong filtering, the vertical jets redistribute the streamwise momentum in the region situated at the bottom of the inertial layer, where these jets reorient themselves to become parallel to the streamwise direction. The end effect is a less smooth increase of $\langle \overline{u}\rangle$ with z compared with cases with small, or no, mussel filtering.

Finally, the random mussel distribution on the channel bed has a negligible effect on the double-averaged streamwise velocity profiles as illustrated in figure 7(e), which compares the predicted profiles from two simulations conducted with BC = 0.056, $ d_{50}/h $ = 0.13 and VR = 0.5 in which a random and, respectively, a staggered regular distribution of the mussels was used. Though significant differences are observed between the local mean vertical velocity profiles, the double-averaged profiles are quite insensitive to the way the mussels are placed on the bed, provided that the distribution of the shells can still be considered uniform at scales larger than 5 times the average distance between the shells.

3.2.1.2. Analytical model describing the vertical distribution of the double-averaged streamwise velocity in fully developed open-channel flow over a mussel bed

Following Yang et al. (Reference Yang, Sadique, Mittal and Meneveau2016) and Wu & Constantinescu (Reference Wu and Constantinescu2025), an analytical model is proposed to approximate the distribution of $\langle \overline{u}\rangle$ inside the channel. The main difference with the model proposed by Wu & Constantinescu (Reference Wu and Constantinescu2025) is that the boundary layer width is replaced by the flow depth, D, in the present work that focuses on the fully developed flow regime. As shown in figure 8(a), the inner, or roughness, layer is divided into a linear sublayer and an exponential sublayer.

Figure 8.

Comparison between predicted vertical profile of $\langle \overline{u}\rangle$ and analytical model. (a) Sketch showing the linear and exponential sublayers of the inner layer and the inertial layer; (b) numerically predicted profiles of $\langle \overline{u}\rangle$ / $ U_{0} $ and analytical model predictions for the BC = 0.280, VR = 0.5, $ d_{50}/h $ = 0.13 simulation; (c) velocity profiles in wall coordinates. Also shown in panels (b) and (c) are the predicted velocity profiles for a channel with uniformly distributed bed roughness ( $ d_{50}/h $ = 0.13) and no mussels.

In all simulations conducted with a mussel bed, the velocity increases linearly with z near the bed inside a region that is much larger than the viscous sublayer developing in cases with $ d_{50}/h $ = 0. Then, the velocity $\langle \overline{u}\rangle$ transitions rapidly to an exponential variation starting at z = δ _m , which identifies the boundary between the two sublayers. The velocity inside the linear sublayer (i.e. $ z \lt \delta_{m}$ ) is approximated as

(3.1) \begin{align}\left\langle \overline{u}\right\rangle '=\frac{U_{m}z}{\delta _{m}},\end{align}

where U _m is the value of $\langle \overline{u}\rangle$ at z = δ _m and the prime symbol indicates the approximate variation assumed by the analytical model. Such a variation is consistent with model 3 of Nikora et al. (Reference Nikora, Koll, McEwan, McLean and Dittrich2004) for the interfacial sublayer that assumes that momentum transfer within the inner (roughness) layer is due to fluid motions comparable with, or even larger than, the height of the roughness elements.

Inside the exponential sublayer ( $ \delta_{m} \lt z \lt h$ ), the velocity is approximated as

(3.2) \begin{align}\left\langle \overline{u}\right\rangle '=U_{h}{\rm e}^{{a(z-h)}/{h}},\end{align}

where U _h is the value of $\langle \overline{u}\rangle$ at the top of the inner layer (z = h given that the protrusion height for all mussels in the array is constant) and a is the so-called attenuation parameter that is a function of the density and spatial distribution of the roughness elements (Yang et al. Reference Yang, Sadique, Mittal and Meneveau2016). An exponential variation is consistent with model 2 of Nikora et al. (Reference Nikora, Koll, McEwan, McLean and Dittrich2004) for the interfacial sublayer that assumes the mixing length within the inner (roughness) layer is constant. Such a model describes the variation of $\langle \overline{u}\rangle$ in deeply submerged plant canopies and, directly relevant for the present study, over the top part of the inner (roughness) layer for cases with tall roughness elements. Using different arguments, the model proposed by Cionco (Reference Cionco1965) also assumes an exponential variation of $\langle \overline{u}\rangle$ inside the inner layer. The fact that, for the particular type of roughness analysed in the present study, a combination of models needs to be used to describe the variation of $\langle \overline{u}\rangle$ over the full height of the inner layer is also consistent with the findings of Nikora et al. (Reference Nikora, Koll, McEwan, McLean and Dittrich2004).

For fully developed flow in an open channel, we assume that the inertial layer extends from z = h to z = D, which formally means that the flow depth D replaces the thickness of the boundary layer induced by the mussels in the fully developed flow region; see also Das, Balachandar & Barron (Reference Das, Balachandar and Barron2022). In this region, the velocity profile is assumed to follow a modified log-law variation supplemented by a law-of-the-wake component (Yang et al. Reference Yang, Sadique, Mittal and Meneveau2016):

(3.3) \begin{align}\left\langle \overline{u}\right\rangle '=\frac{u_{\tau }}{\mathit{\kappa }}\left[ln\left(\frac{z-d_{0}}{z_{0}}\right)+ \Pi w\left(\frac{z}{D}\right)\right].\end{align}

Here u _τ is the bed friction velocity, κ is the von Karman constant, d ₀ is the displacement height (Jackson Reference Jackson1981), z ₀ is the hydrodynamics roughness length and Π is the scaling coefficient of the wake function. Following Coles (Reference Coles1956), the wake function used was w( $ z/D = 2\text{sin}^{2}[\pi z/(2D)] $ . In the present simulations, d ₀ is estimated as the distance from the bed where the resultant drag force is applied on the protruding part of each mussel (Wu & Constantinescu Reference Wu and Constantinescu2025). Basically, (3.3) modifies the origin of z in the standard logarithmic law used for rough surfaces by d ₀. The values of u _τ and z ₀ are obtained by imposing $\langle \overline{u}\rangle '=\langle \overline{u}\rangle$ at the bottom and top of the inertial layer. The scaling coefficient is used as a tuning parameter to obtain the overall best agreement between $\langle \overline{u}\rangle '$ and $\langle \overline{u}\rangle$ .

Simulations performed with a mussel bed (BC = 0.280) and with no mussels (BC = 0) and distributed bed roughness ( $ d_{50}/h $ = 0.13) are used to illustrate how the presence of sparse, large roughness elements modifies the vertical profile of $\langle \overline{u}\rangle$ and how well the analytical model given by (3.1)–(3.3) can approximate the profile determined based on 3-D simulations. As shown in figure 8(b), the presence of the mussels and the associated reduction of $\langle \overline{u}\rangle$ inside the inner layer is accompanied by a change in the shape of the profile inside this layer. Equations (3.1) and (3.2) provide a good approximation of the predicted velocity profile inside the inner layer, while the boundary between the linear and exponential sublayers occurs around $ h/3$ (figure 8 b). An important observation is that the law-of-the-wall component does not provide a good approximation of the velocity profile over most of the inertial layer (e.g. large differences are observed between the predicted profile of $\langle \overline{u}\rangle$ and the law-of-the-wall component for $ z/D \gt 0.3 $ in figure 8 b). However, when supplemented by a law-of-the-wake component with a properly fitted value of Π, the agreement between the predicted profile of $\langle \overline{u}\rangle$ and the analytical model given by (3.3) is very good. Another important observation is that, when plotted in wall coordinates, the double-averaged velocity profile for the simulation conducted with no mussels shows the correct behaviour, with a log-law region that is parallel to that expected for a smooth wall but with values of U ⁺ = $\langle \overline{u}\rangle/u_{\tau}$ that are shifted downwards (i.e. ΔU ⁺ ≈ 9.0 for $ z^{+} = zu_{\tau}/\nu \gt 200$ in figure 8(c), where ν is the molecular viscosity). When protruding mussels are present at the bed surface, the profile of U ⁺ does not display the usual log-law region mostly because of the large contribution of the law-of-the-wake component.

3.2.1.3. Effects of varying BC, d₅₀/h and VR on the main parameters of the analytical model describing the streamwise velocity variation, the bed friction velocity and the equivalent bed roughness height

Table 1 summarises the main parameters and variables of the analytical model used to approximate $\langle \overline{u}\rangle .$ Though both U _m and δ _m are generally affected by the values of BC and $ d_{50}/h $ (table 1), the non-dimensional slope of the velocity gradient inside the linear sublayer, $ (U_{m}/U_{0})/(\delta_{m}/h)$ , varies in a narrow range (i.e. 0.13–0.18) for BC ≥ 0.056. As the siphons are situated at the top of the inner layer, the influence of the velocity ratio on the velocity profile inside the linear sublayer is negligible. For constant $ d_{50}/h $ , the thickness of the linear sublayer decreases with increasing BC, which also means that the thickness of the exponential sublayer increases with BC. This result is consistent with the trend reported by Wu & Constantinescu (Reference Wu and Constantinescu2025) for developing, mussel-induced, boundary layers in an open channel, as well as with findings from several experimental studies (Nikora et al. Reference Nikora, Green, Thrush, Hume and Goring2002; Quinn & Ackerman Reference Quinn and Ackerman2014). However, in the present simulations of fully developed flow, the rate of decrease of δ _m is very small for 0.056 ≤ BC ≤ 0.280, which is the range of interest for most natural streams containing mussels. Meanwhile, the effect of increasing the distributed bed roughness for a fixed value of BC is to increase the thickness of the linear sublayer. The increase of the thickness of the exponential sublayer with BC is explained by the fact that a larger part of the bottom of the mussels are situated in the wake generated by upstream mussels as the distance between the mussels decreases. In other words, the sheltering height increases with BC. For sufficiently high surface mussel coverage densities (i.e. for BC = 0.8), the sheltering height approaches the height of the roughness element (i.e. h) and transition to the skimming flow regime is complete.

As BC increases, the decay of $\langle \overline{u}\rangle$ inside the exponential sublayer will be larger with respect to the velocity at the top of the roughness layer, U _h . This is the main reason for the observed increase of the attenuation parameter a with BC for constant $ d_{50}/h $ and VR (see figure 9(a) and table 1). The increase of a with BC for constant $ d_{50}/h $ is accompanied by a decrease of U _h (table 1). Meanwhile, U _h is fairly insensitive to variations in VR and $ d_{50}/h $ (table 1). The increase of δ _m with $ d_{50}/h $ for constant BC is the main reason for the observed monotonic increase of a with $ d_{50}/h $ in figure 9(b). Not surprisingly, the effect of the velocity ratio, VR, is negligible on the flow inside the exponential sublayer (figure 9 c). It is also relevant to mention that Yang et al. (Reference Yang, Sadique, Mittal and Meneveau2016) predicted a minimum value a = 0.4 for a boundary layer over cubic roughness elements for cases with negligible sheltering. The minimum value in the simulation with the lowest value of the surface mussel coverage density (i.e. BC = 0.015) is 0.5 (table 1). In the case of a mussel-induced boundary layer developing in an open channel, Wu & Constantinescu (Reference Wu and Constantinescu2025) estimated significantly lower values of a for comparable BC. This is likely due to the lower values of $ U_{h}/U_{0}$ predicted for cases when such boundary layers develop in a depth-limited flow where the flow outside the boundary layer accelerates in the streamwise direction.

Figure 9.

Attenuation factor a (left) and equivalent (sand-grain) roughness height in wall units $ K_{S} $ ⁺ (right) as a function of (a) BC in the simulations conducted with $ d_{50}/h $ = 0.13 and VR = 0.5, (b) $ d_{50}/h $ in the simulations conducted with BC = 0.112 and VR = 0.5, (c) VR in the simulations conducted with BC = 0.112 and $ d_{50}/h $ = 0.13.

The values of the wake parameter in the present simulations conducted with constant h are mainly a function of BC (table 1). Until the transition to the skimming flow regime starts, Π is increasing monotonically with BC, as the increase in the surface mussel coverage density accelerates the flow near the top (free surface) boundary and decelerates the flow near the bottom of the inertial layer (figure 7 a). Though this finding is consistent with the trend observed for a developing mussel-induced boundary layer in an open channel, the values of Π are much smaller in the case of fully developed flow (i.e. 0.15 ≤ Π ≤ 0.3) compared with those (0.55 ≤ Π ≤ 1.25) reported by Wu & Constantinescu (Reference Wu and Constantinescu2025) for arrays with a similar level of mussel burrowing. The large decrease of Π between developing and fully developed open-channel flows (e.g. from 1.25 to 0.2 for BC = 0.056) is mainly because in a developing boundary layer in a channel the flow at the edge of the boundary layer is driven by the accelerating flow in the ‘free-stream’ region, which induces a large wake component. This effect is not present in the case of a fully developed flow, which also explains why the range of Π values in the present simulation is close to typical values (Π ≈ 0.2) for boundary layers over large-scale roughness elements where the free-stream velocity is close to constant (Yang et al. Reference Yang, Sadique, Mittal and Meneveau2016) and in turbulent open-channel flows (e.g. Nezu & Rodi Reference Nezu and Rodi1986; Kirkgöz Reference Kirkgöz1989). Once the flow over the mussel bed transitions to the skimming flow regime, Π decreases as the mussel bed starts behaving more like uniformly distributed roughness.

The double averaging formally removes the mussels and, if present, the distributed roughness from the computational domain. However, the bed friction velocity, u _τ , calculated based on the double-averaged velocity profile, $\langle \overline{u}\rangle$ (z), accounts for the presence of these two types of roughness. As expected, for constant surface mussel coverage density, $ u_{\tau}/U_{0}$ , increases with $ d_{50}/h $ and is not affected by the velocity ratio, VR. Until the skimming flow regime is reached, $ u_{\tau}/U_{0} $ increases slightly with BC (table 1), as the growth in the number of large-scale roughness elements per unit area translates into an increase of the equivalent bed roughness, as will be discussed later. For all mussel bed simulations, the ratio $ u_{\tau}/U_{D}$ varies little with both BC and $ d_{50}/h $ (table 1).

Figure 10.

Non-dimensional vertical profiles of $\langle \overline{u}\rangle$ inside the inertial layer in wall coordinates (log-linear scale) used to estimate the values of the equivalent (sand-grain) roughness height, $ K_{S} $ , of the log-law component of $\langle \overline{u}\rangle$ for z > > d ₀. (a) Effect of varying BC in the simulations performed with $ d_{50}/h $ = 0.13 and VR = 0.5; (b) effect of varying $ d_{50}/h $ in the VR = 0.5 simulations performed with BC = 0.056 and BC = 0.280. The dashed–dotted lines show the standard law of the wall in a fully developed open-channel flow over a rough bed with sand-grain roughness $ K_{S} $ .

Following Wu & Constantinescu (Reference Wu and Constantinescu2025), the profiles of $\langle \overline{u}\rangle$ are used to determine the equivalent (sand-grain) roughness height, $ K_{S} $ . In the absence of the law-of-the-wake component and assuming the logarithmic layer is situated sufficiently far from the channel bed, i.e. z > > d ₀, (3.3) reduces to the standard form of the law of the wall for a rough surface, i.e.

(3.4) \begin{align}\left\langle \overline{u}\right\rangle =u_{\tau }\left[\frac{1}{\mathit{\kappa }}ln\left(\frac{z}{K_{S}}\right)+B_{S}\right],\end{align}

where B _S is the law-of-the-wall constant. Figure 10 includes the law of the wall and the associated value of $ K_{S} $ ⁺ = $ K_{S} $ $ u_{\tau}/U_{0}$ for fully developed flow over a rough bed with uniformly sized, sand-grain roughness. For the rough-bed surface with $ d_{50}/h $ = 0.13, the inferred value of $ K_{S} $ ⁺ is 176. This corresponds to $ K_{S}/d_{90} $ = 2.6, which is close to standard values of this ratio assumed in rivers (i.e. $ K_{S}/d_{90} $ = 2.0–3.0). As expected, increasing the roughness of the bed on which the mussels are placed results in an increase of $ K_{S} $ . Results in figures 9(b) and 10(b) show that this effect decreases at high surface mussel coverage densities (e.g. the increase of $ K_{S} $ ⁺ is of the order of 400 % for BC = 0.056 but only of about 40 % for BC = 0.280 if one compares simulations conducted with a smooth bed surface and with $ d_{50}/h $ = 0.13). This effect is expected to be negligible once the skimming flow regime is observed as the flow inside the inertial layer is not anymore affected by the smaller-scale distributed bed roughness. For comparable conditions, $ K_{S} $ ⁺ has similar values for fully developed channel flow and for a developing boundary layer in a channel. A direct comparison between present results and the cases investigated by Wu & Constantinescu (Reference Wu and Constantinescu2025) can be made only for a mussel bed with BC = 0.056, $ h/H $ ≈ 0.5 and $ d_{50}/h $ = 0 for which $ K_{S} $ ⁺ = 76 for fully developed flow and $ K_{S} $ ⁺ = 60 for a developing boundary layer. These values correspond to $ K_{S}/h $ ≈ 1. If mussels are placed on a rough-bed surface, $ K_{S}/d_{90} $ increases with respect to the limiting case where no mussels are present due to the additional ‘roughness’ provided by the mussels. For example, for $ d_{50}/h $ = 0.13, the predicted $ K_{S}/d_{90} $ ranges between 5 and 8 for 0.056 ≤ BC ≤ 0.280 (table 1), which is about 2 times the value ( $ K_{S}/d_{90} $ ≈ 3) predicted in the corresponding simulation with no mussels.

For constant BC and $ d_{50}/h $ , $ K_{S} $ ⁺ decreases with increasing velocity ratio (figure 9(c) and table 1), though the decrease is only of the order of 10 % for 0 ≤ VR ≤ 1.22 (table 1). This effect is basically controlled by the vertical momentum of the excurrent siphon jets.

For constant $ d_{50}/h $ and VR, $ K_{S} $ ⁺ increases with BC until the skimming flow regime is reached (see table 1 and figure 10 a), a result similar to that observed for developing boundary layers in open channels with a mussel bed (Wu & Constantinescu Reference Wu and Constantinescu2025). This is mainly due to the increase of the average drag force per streamwise length acting on the fluid inside the inner layer. As the flow transitions to the skimming flow regime, $ K_{S} $ ⁺ peaks before starting to decrease with increasing BC. This is because as mussels get closer to each other, the bottom part of the inner layer contributes negligibly to the drag force acting on the mussels, so the ‘effective height’ of the large-scale roughness elements starts decreasing with BC. Assuming a smooth bed, $ K_{S} $ ⁺ in the BC = 0.828 simulation is expected to be comparable with the value predicted for a mussel bed with BC ≈ 0.224 and $ d_{50}/h $ = 0. Another interesting question is what should be the size of the sediment (e.g. d ₉₀) for the gravel bed to match the $ K_{S} $ for a mussel bed with a given BC and $ d_{50}/h $ = 0. For example, present results show that a mussel bed with BC ≈ 0.15 containing mussels with a relative protrusion height $ h/H $ = 0.48 placed on a smooth bed will have the same $ K_{S} $ (≈ 2.6d ₉₀) as a gravel bed with $ d_{50}/h $ = 0.13 ( $ d_{90}/h$ = 0.21).

3.2.2.1. Effects of varying BC, d₅₀/h, VR and the mussels’ arrangement on the bed surface

A general characteristic of the vertical profile of the double-averaged TKE in the simulations conducted with 0.056 ≤ BC ≤ 0.280 is the presence of two relative maxima, one inside the inner layer and one at the bottom of the inertial layer (figure 11 a). The peak inside the inner layer is likely associated with the generation of turbulent eddies inside the horizontal separated shear layers on the two sides of the protruding mussel and vortex shedding. The peak inside the inertial layer is primarily due to eddies shed inside the vertical separated shear layer near the top of the protruding mussels. By contrast, only one peak situated just above the inner layer is observed in mussel bed simulations where the skimming flow regime dominates (see, e.g. the results for BC = 0.828 in figure 11 b) or in simulations with distributed bed roughness but no mussels (see, e.g. the results for BC = 0 in figure 11 a). This is because for sufficiently high mussel clustering, vortex shedding and generation of energetic eddies inside the horizontal separated shear layers is suppressed. However, energetic eddies are still generated in the vertical separated shear layers (see, e.g. figure 5 c). The largest $\langle \textit{TKE} \rangle$ corresponds to the inner layer peak for low BC values. For BC ≥ 0.112, the maximum $\langle \textit{TKE} \rangle$ occurs at the bottom of the inertial layer. Overall, the largest levels of amplification of $\langle \textit{TKE} \rangle$ inside the water column are generated in the simulations with the largest values of $ K_{S} $ ⁺ . Consistent with this, $\langle \textit{TKE} \rangle$ decreases once the skimming flow regime is established (figure 11 b).

Figure 11.

Vertical profiles of $\langle \textit{TKE} \rangle$ as a function of (a) BC in the simulations conducted with $ d_{50}/h $ = 0.13 and VR = 0.5; (b) BC in the simulations conducted with $ d_{50}/h $ = 0.0; (c) $ d_{50}/h $ in the simulations conducted with VR = 0.5, BC = 0.056 and BC = 0.280; (d) VR in the simulations conducted with BC = 0.280 and $ d_{50}/h $ = 0.13; (e) random versus irregular distribution of the mussels in the simulations conducted with BC = 0.056, $ d_{50}/h $ = 0.13 and VR = 0.5. Frame (f) compares the resolved and modelled components of the $\langle \textit{TKE} \rangle$ in the simulation conducted with BC = 0.056, $ d_{50}/h $ = 0.00 and VR = 0.5. The black and red arrows in the frame point toward the $\langle \textit{TKE} \rangle$ peaks inside the inner and inertial layers, respectively.

The presence of distributed bed roughness in the mussel bed simulations results in an increase of $\langle \textit{TKE} \rangle$ over the whole water column mainly due to the generation of energetic eddies by the gravels. However, this effect becomes negligible for BC ≥ 0.280 (figure 11 c). The main effect of increasing the mussel filtering discharge is an increase of $\langle \textit{TKE} \rangle$ inside the inertial layer (figure 11 d).

The use of a regular staggered distribution of the mussels has a negligible effect on the vertical distribution of $\langle \textit{TKE} \rangle$ compared with the case when a randomised distribution is used (see, e.g. figure 11(e) that compares results for the BC = 0.056, $ d_{50}/h $ = 0.13, VR = 0.5 case). In particular, the peak $\langle \textit{TKE} \rangle$ values inside the inner and inertial layers are very close.

Given the use of DES, it is also relevant to check the relative contribution of the modelled part of $\langle \textit{TKE} \rangle$ . As an example, figure 11(f) compares the modelled and resolved parts of $\langle \textit{TKE} \rangle$ for the BC = 0.056, $ d_{50}/h $ = 0, VR = 0.5 case. The modelled part is close to zero inside the inertial layer and over most of the exponential sublayer. The contribution of the modelled part is significant only very close to the bed surface, but the ratio between the peak values of the resolved and model parts is larger than 20. So, in a good approximation, the resolved part of the turbulence kinetic energy approximates the total value of $\langle \textit{TKE} \rangle$ . This is why the modelled part of $\langle \textit{TKE} \rangle$ was not plotted in figure 11(a–e).

3.2.2.2. Analytical model describing the vertical distribution of $\langle \textit{TKE} \rangle$ inside the inertial layer for fully developed open-channel flow over a mussel bed

In the case of a mussel-induced developing boundary layer in a smooth bed open channel, Wu & Constantinescu (Reference Wu and Constantinescu2025) have shown that the scaled $\langle \textit{TKE} \rangle$ profiles inside the inertial layer are close to self-similar above the vertical location, $z_{\textit{TKE},\textit{max}}$ , where the maximum $\langle \textit{TKE} \rangle$ inside the inertial layer, TKE _max , is observed. Inside most of the inertial layer, the decay of the $\langle \textit{TKE} \rangle$ was found to be close to exponential. As the top of the boundary layer was approached, the decay was better approximated by a linear function. The values of $\langle \textit{TKE} \rangle / \textit{TKE}_{\textit{max}}$ at the edge of the boundary layer were of the order of 0.02–0.04. The outer turbulent flow corresponded to the fully developed open-channel flow over a smooth bed. Present results obtained for a fully developed flow over a mussel bed also suggest an exponential decay of the $\langle \textit{TKE} \rangle$ , as exemplified by results in figure 12(a) for a simulation conducted with BC = 0.112, $ d_{50}/h $ = 0.13 and VR = 0.5. In the case of a fully developed flow the boundary layer width is replaced by the channel depth, D, and the scaling of the vertical coordinate is defined as

Figure 12.

Scaled profiles of $\langle \textit{TKE} \rangle$ /TKE _max inside the inertial layer (z > z _TKE ): (a) BC = 0.112, $ d_{50}/h $ = 0.13, VR = 0.5; (b) BC = 0, $ d_{50}/h $ = 0.13, VR = 0.5; (c) effect of varying BC for cases with $ d_{50}/h $ = 0.13 and VR = 0.5. The dashed lines show the DES predictions while the solid lines show the analytical model predictions.

(3.5) \begin{align}{z}_{\textit{TKE}}^{\prime}=\frac{z-z_{\textit{TKE},\textit{max}}}{D-z_{\textit{TKE},\textit{max}}}.\end{align}

As opposed to the developing boundary layer cases, the $\langle \textit{TKE} \rangle$ decay is close to exponential until the free surface and the values of $\langle \textit{TKE} \rangle / \textit{TKE}_{\textit{max}} $ at z = D are about one order of magnitude higher (figure 12 c). The absence of the linear decay region and the higher values of $\langle \textit{TKE} /TKE_{\textit{max}} \rangle$ are easily explained by the fact that the flow inside the boundary layer does not have to transition to a low-turbulence level outer flow. Moreover, the same exponential decay of $\langle \textit{TKE} \rangle / TKE_{\textit{max}} $ applies for cases with a gravel bed but with no mussels, as exemplified in figure 12(b) for a simulation conducted with $ d_{50}/h $ = 0.13 and VR = 0.5.

To better examine trends on how the scaled $\langle \textit{TKE} \rangle$ varies inside the inertial layer with the main parameters describing the bed surface, one can compare the profiles given by the analytical model approximating the decay of $\langle \textit{TKE} \rangle / \textit{TKE}_{\textit{max}} $ :

(3.6) \begin{align}\frac{\langle \textit{TKE} \rangle }{\textit{TKE}_{max}}={\rm e}^{{r_{\textit{TKE}}}{z}_{\textit{TKE}}^{\prime}}.\end{align}

Here r _TKE is the parameter of the exponential decay determined such that the best overall agreement is observed with the profile inferred from the corresponding simulation. Figure 12(c) shows the effect of varying BC in the simulations conducted with a gravel bed ( $ d_{50}/h $ = 0.13). The value of $\langle \textit{TKE} \rangle / \textit{TKE}_{\textit{max}} $ at the free surface increases monotonically with BC. This increase is mostly related to the increase of the number of energetic eddies shed by the top parts of the mussels. Compared with the limiting case of BC = 0, $\langle \textit{TKE} \rangle /\textit{TKE}_{\textit{max}} $ at the free surface increases by about 40 % in the BC = 0.280 simulation (figure 12 c). The increase basically stops once the skimming flow regime becomes dominant. The decays of $\langle \textit{TKE} \rangle / \textit{TKE}_{\textit{max}}$ with increasing $z'_{\textit{TKE}}$ in the simulations conducted with BC = 0.280 and BC = 0.828 are close to identical.

Figure 13 summarises the main trends in the variations of the peak $\langle \textit{TKE} \rangle$ inside the inertial layer, its vertical position, $ z_{\textit{TKE},\textit{MAX}}$ , and of the exponential decay parameter in the analytical model for cases where mussels can still be characterised as sparse roughness elements (i.e. BC ≤ 0.280). The increase of $ z_{\textit{TKE},\textit{MAX}}/h$ with BC is close to linear and independent of the presence of smaller-scale distributed roughness at the channel bed. For both a smooth bed and a gravel bed, $ \textit{TKE}_{\textit{max}} / U_{0}^{2}$ increases with BC but the rate of increase is larger for the smooth bed cases such that the values are about the same at values of BC where the proximity of the mussels is sufficiently high for the turbulence generated by the gravels to have a negligible effect on the TKE inside the inertial layer (i.e. BC = 0.280 in the present simulations). The main result of increasing BC is to increase the rate of decay of $\langle \textit{TKE} \rangle / \textit{TKE}_{\textit{max}} $ (i.e. r _TKE ) in the presence of a rough bed (figure 13 c). The effect is opposite in the case of a smooth bed. As for $ \textit{TKE}_{\textit{max}} / U_{0}^{2}$ , the effect of the presence of smaller-scale distributed roughness on r _TKE becomes negligible for BC = 0.280.

Figure 13.

Effect of varying BC on the non-dimensional variables and parameters in the analytical model approximating the vertical variation of the $\langle \textit{TKE} \rangle$ inside the inertial layer: (a) $ z_{\textit{TKE},\textit{MAX}}$ , (b) TKE _max and (c) r _TKE . Results are compared for simulations conducted with $ d_{50}/h $ = 0 and 0.13 and VR = 0.5.

3.2.3.1. General features of the shear stress distributions in fully developed channel flow over a mussel bed

Several studies of open-channel flows and developing boundary layers over a rough bed (e.g. Forooghi et al. Reference Forooghi, Stroh, Schlatter and Frohnapfel2018; Jelly & Busse Reference Jelly and Busse2019; Wu & Constantinescu Reference Wu and Constantinescu2025) have shown that dispersive stresses can have a significant contribution to the total stress inside the inner layer and its immediate vicinity. The vertical profiles of the primary Reynolds shear stress, $\langle \overline{u^{\prime}w^{\prime}}\rangle$ , and of the corresponding dispersive stress, $\langle (\overline{u}-\langle \overline{u}\rangle )(\overline{w}-\langle \overline{w}\rangle )\rangle$ , are compared in figure 14 to illustrate the effect of varying BC and $ d_{50}/h $ in the present mussel bed simulations.

Figure 14.

Vertical profiles of the primary resolved Reynolds (solid lines) and dispersive (dashed) shear stresses. (a) Effect of varying BC in the simulations conducted with $ d_{50}/h $ = 0 and VR = 0.5; (b) effect of varying $ d_{50}/h $ in the simulations conducted with BC = 0.056 and VR = 0.5; (c) random versus irregular distribution of the mussels in the simulations conducted with BC = 0.056, $ d_{50}/h $ = 0.13 and VR = 0.5. Frame (d) compares the resolved (solid lines) and modelled (dash-dotted lines) components of the double-averaged primary shear stress in the simulations conducted with BC = 0.056, VR = 0.5 $ d_{50}/h=0 $ and $ d_{50}/h =0.13$ .

A general feature of the profiles of the resolved shear Reynolds stress $\langle \overline{u^{\prime}w^{\prime}}\rangle$ in figure 14 is the presence of two peaks, the main one situated at the bottom of the inertial layer and the secondary one situated inside the inner layer, around $ z/h $ = 0.7. Such a double-peak distribution was not present in the developing boundary layer simulations of Wu & Constantinescu (Reference Wu and Constantinescu2025). However, there is an important difference between the two series of simulations. In the present simulations the main axis of the mussels makes an angle of about $45^\circ$ with the vertical direction as the mussels are tilted upstream. By contrast, the main axis of the mussels was close to vertical in the simulations conducted by Wu & Constantinescu (Reference Wu and Constantinescu2025) such that the asymmetry between the front and back sides of the protruding mussel was fairly low. Moreover, both the fully developed channel flow simulations of Jelly & Busse (Reference Jelly and Busse2019) conducted with irregular near-Gaussian roughness elements and the present simulations conducted with no mussels and distributed irregular roughness with $ d_{50}/h \gt 0 $ predicted that the vertical profiles of $\langle \overline{u^{\prime}w^{\prime}}\rangle$ contain only one peak situated close to the interface between the inner and the inertial layers. Finally, as for the TKE, the contribution of the modelled shear stress to the total value of the double-averaged primary shear stress is negligible except very close to the bed surface. This result is illustrated in figure 14(d) for the smooth bed and rough-bed simulations conducted with BC = 0.056 and VR = 0.5.

The dispersion shear stress profiles predicted in the present simulations with BC ≥ 0.056 also show one qualitative difference with those observed for a developing boundary layer in an open channel containing a mussel bed. The peak of the dispersion shear stress in figure 14 always occurs at the interface between the inner layer and the inertial layer ( $ z/h $ ≈ 1). By contrast, the corresponding simulations of Wu & Constantinescu (Reference Wu and Constantinescu2025) predicted the peak to occur inside the inner layer ( $ z/h $ = 0.6–0.7). This difference is not necessarily due to the different type of flow over the mussel bed (i.e. fully developed versus developing) as the peak of the dispersion stress was situated at $ z/h $ = 0.5–0.65 in channel flow simulations of Jelly & Busse (Reference Jelly and Busse2019) conducted with irregular, closely spaced, near-Gaussian roughness elements. Rather, this qualitative difference is most likely due to the large inclination of the partially burrowed mussels in the present simulations. This is also supported by the fact that qualitative changes in the dispersion stress profiles were observed in fully developed channel flow over other types of roughness characterised by large irregularities in terms of the shape and/or local degree of clustering of the roughness elements. For example, the dispersion shear stress profiles predicted by Sarakinos & Busse (Reference Sarakinos and Busse2024) for fully developed open-channel flow over a smooth surface containing barnacle clusters covering 10 % of the bed surface showed the dispersion stress contained a negative peak above the vertical elevation where the main positive peak occurred.

3.2.3.2. Effects of varying BC, d₅₀/h and the mussels’ arrangement on the bed surface

As exemplified in figure 14(a), for constant $ d_{50}/h $ , both the Reynolds and the dispersive shear stresses are amplified with increasing surface mussel coverage density. The ratio between the peak values of the Reynolds and dispersive stresses is fairly insensitive to the increase in BC. By contrast, increasing the roughness of the bed surface on which the mussels are placed results in an increase of the Reynolds shear stress over the entire water column (e.g. for BC = 0.056, the increase is of the order of 50 % for the $ d_{50}/h $ = 0.13 case compared with the $ d_{50}/h $ = 0 case), but has basically no effect on the dispersive stress (figure 14 b). This is because the dispersive stresses are mostly associated with flow non-uniformity and secondary flow generated by the larger-scale roughness elements (i.e. protruding mussels in the present simulations) present at the bed surface. Meanwhile, the distributed bed roughness generates lots of energetic turbulent eddies (see, e.g. figures 5 a and 5 b) that increase both the normal and shear Reynolds stresses.

It is also relevant to compare the peak values of the non-dimensional Reynolds and dispersive shear stresses for fully developed and developing open-channel flow simulations conducted with the same value of BC and with mussels placed on a smooth surface ( $ d_{50}/h $ = 0). Such a comparison is done for BC = 0.056 and $ h/H $ ≈ 0.5. For a developing boundary layer in an open channel, Wu & Constantinescu (Reference Wu and Constantinescu2025) reported peak absolute values of the dispersive and shear stresses of 0.004 ${U}_{0}^{2}$ and 0.008 ${U}_{0}^{2}$ , respectively, at distances where the thickness of the boundary layer was close to the channel depth. The corresponding peak values in the present simulations are 0.001 ${U}_{0}^{2}$ and 0.0045 ${U}_{0}^{2}$ , respectively (figure 14 a). So, the transition to the fully developed regime is accompanied by a relative reduction of the Reynolds and dispersive stresses and by an increase of the ratio of the peak value of the Reynolds and dispersive stresses. For BC = 0.056 and $ h/H $ ≈ 0.5, this ratio is about 2.0 for the developing boundary layer flow and close to 4.5 for the fully developed flow with $ K_{S} $ ⁺ = 76. The latter result is close to the value of about 4.0 observed in fully developed flow over irregular, near-Gaussian roughness elements with comparable values of $ K_{S} $ ⁺ by Jelly & Busse (Reference Jelly and Busse2019).

The effect of varying the spatial distribution of the mussels on the double-averaged profiles of the primary Reynolds and dispersive shear stresses is negligible (see, e.g. figure 14(c) for the simulations conducted with BC = 0.056, $ d_{50}/h $ = 0.13 and VR = 0.5).