Stabilization of fluvial bed sediments by invasive quagga mussels (Dreissena bugensis)

River gravel mobility is an important control on river behaviour, morphology, and ecosystem processes. Gravel stability is dependent on abiotic flow and sediment properties, alongside less widely acknowledged biotic processes. The quagga mussel (Dreissena bugensis), a highly invasive bivalve, frequently occurs at high population densities in rivers and lakes. Quagga mussels attach to benthic sediments using byssal threads, which might affect sediment stability and thereby broader river geomorphology. We aimed to (1) characterize controls of quagga mussel sediment attachment by conducting a field survey in an invaded river, (2) investigate resultant changes in the critical shear stress needed to entrain fluvial bed materials via an ex situ flume experiment, for measured average (135 m−2) and potential future (270 m−2) mussel densities, and (3) model how this may affect sediment transport rates. From field surveys, mean quagga mussel density was 122 m−2, attaching to an average of 591 g m−2 of mineral bed sediments. Across the survey reach, mussels attached to all grain sizes available, with attachment driven by grain availability, rather than active selection of particular grain sizes. In the ex situ flume experiment, densities of 135 mussels m−2 did not significantly increase the critical shear stress of fluvial bed materials, but a density of 270 mussels m−2 significantly increased critical shear stress by 40%. Estimates of the proportion of time these critical stresses are exceeded at the field site indicated high densities of quagga mussels may reduce the occurrence of a geomorphically active flood event from Q30 to Q2. These results present feasible invasion scenarios, as quagga mussels frequently reach benthic densities orders of magnitude greater than observed here. This study indicates that substantial alterations to bedload sediment transport may occur following quagga mussel invasion. Future geomorphic modelling should include biology to better understand rates and processes of landscape development.

Such stabilization may have important geomorphic effects. For example, Cardinale et al. (2004)  Whilst many species alter fluvial sediment regimes, invasive species are worthy of particular attention. Population density and exhibited behavioural traits are key controls on the zoogeomorphic potential of a species (Albertson & Allen, 2015;Mason & Sanders, 2021). Invasive species frequently achieve very high population densities in invaded habitats (Hansen et al., 2013) and may develop novel zoogeomorphic behaviours upon invasion (Reader & Laland, 2003;Sol & Weis, 2019). Invasive species may therefore cause stronger geomorphic disturbances to river systems when compared to native species (sensu Butler, 2006;Crooks, 2002;Emery-Butcher et al., 2020;Fei et al., 2014;Harvey et al., 2011).
Quagga mussels (Dreissena bugensis) are a highly invasive bivalve, from the Ponto-Caspian region, rapidly spreading throughout freshwater systems in Europe and North America, particularly lowland river catchments (e.g., Moselle, Bij de Vaate & Beisel, 2011;Rhine, Imo et al., 2010;Thames, Himson et al., 2020;Upper Mississippi, Grigorovich et al., 2008;Volga, Orlova et al., 2004). As adults, quagga mussels are usually sessile, up to 45 mm long and frequently occur at population densities exceeding 1000 individuals m À2 (Nalepa et al., 2009;Strayer et al., 1999). Notably, quagga mussels secrete byssal threads to attach themselves to the benthic substrate, forming sediment-mussel complexes called 'druses' (Figure 1). The presence of quagga mussel druses in a river have been hypothesized to alter sediment transport through multiple mechanisms (Mills, 2019;Mason & Sanders, 2021); (1) by compressing the underlying riverbed through the mass of the druse, (2) by altering near-bed hydraulics as mussel shells may protrude above the bed, and (3) by attaching particles together, increasing their effective size, and binding them to the riverbed. At high densities, many druses may therefore stabilize the gravel bed by forming large-scale protective 'matting' over the substrate. To our knowledge, only one study has investigated the geomorphic effects of riverbed stabilization by quagga mussels. Mills (2019) used flume experiments to investigate whether the presence of quagga mussels reduced sediment transport when exposed to a flow velocity exceeding the critical entrainment of river gravels. A density of 250 mussels m À2 reduced sediment transport by 59% at a high flow velocity above the critical entrainment stress compared to when mussels were absent, giving the first evidence that quagga mussels may reduce sediment transport. However, this only considered one flow velocity, and so changes in the critical shear stress required to entrain gravels, changes to bedload sediment transport across a flood event, and the frequency of which this may be geomorphically important in the field are unknown. Further, no studies have characterized the geomorphic characteristics of quagga mussel druses in the field. Therefore, the geomorphic characteristics and impacts of byssal thread anchoring by quagga mussels are unknown.
In this study, first, we use a field study of an invaded gravel bed river to characterize the sediment attachment of quagga mussel druses. Second, we measure the change in critical shear stress required to mobilize riverbed gravels in the presence of quagga mussels using flume experiments. Third, we estimate how these changes in critical shear stress might affect the proportion of time that bed materials in the study reach are mobile, and estimate bedload sediment transport rates using the Meyer-Peter and Müller (1948) bedload equation under current and projected invasion scenarios.
Specifically, the following questions were investigated: 1. What are the characteristics of quagga mussel sediment attachment, and what are the controls on sediment attachment in quagga mussel druses? F I G U R E 1 Quagga mussels in the Wraysbury River, UK. Images show individual mussel-sediment complexes or 'druses' (a, b), and an in situ population (c; density: 204 quagga mussel individuals m À2 ). Scale approximate [Color figure can be viewed at wileyonlinelibrary. com] 2. Do quagga mussel druses increase the shear stress for incipient motion and subsequent bedload transport rates of riverbed gravels?
3. Which grain size fractions do quagga mussel druses most effectively stabilize? 4. Do quagga mussels increase the return period for geomorphically active floods? 5. Do quagga mussels reduce bedload sediment transport rates over annual timescales?

| METHODS
We coupled an in situ survey of quagga mussels and druse sediment in an invaded UK river (Q1) with hydraulic flume experiments to examine the mobility of bed sediment under different densities of quagga mussels (Q2 and Q3). Results from both the field study and flume experiments were then used to model bedload sediment transport under different quagga mussel invasion scenarios in the river (Q4), increasing the spatial and temporal scale of the study.

| Characterizing sediment attachment by quagga mussels in the field
The Wraysbury River is an 8.7 km side channel of the River Colne near Staines-upon-Thames, UK (51.4514, À0.5209; Figure 2), which for historical reasons has been separately named. The Wraysbury River was selected as a study site as it is currently the only known wadable river in the UK that is heavily colonized by quagga mussels. Whilst quagga mussels typically prefer deeper habitats (Jones & Ricciardi, 2005;Mills et al., 1996), the distribution of Dreissenid mussels is also dependent on dispersal related factors (Mills, 2019;Quinn et al., 2014;Thorp et al., 2002), and the presence of hard mineral substrates (Karatayev et al. 1998;Strayer et al., 1996). As such, their high density in the Wraysbury River (54 mussels m À2 from 2015 to 2016; Mills et al., 2017) likely reflects the abundance of preferable mineral substrate, and the direct connection with the Wraysbury Reservoir, considered to be the original source of quagga mussels in the UK (Blackman et al., 2020). The Wraysbury River therefore represents an appropriate site to understand the potential impacts of quagga mussels on fluvial processes on other similar rivers given their anticipated widespread distribution throughout the River Thames catchment, throughout the UK, and throughout Europe and North America (e.g., Gallardo & Aldridge, 2018;Quinn et al., 2014).
The morphology of the Wraysbury River is homogenous, straightened and over-deepened with poorly developed bed topography throughout its length, with a consistent width of 5 to 7 m and depth of < 0.5 m. Mean river width at the study site was 6.4 m, mean depth was 0.35 m, and mean bankfull depth was 0.97 m. Measured morphological features are representative of the rest of the river channel adjacent to Staines Moor. Channel thalweg slope was measured by levelling to be 0.47% for the length of the 60 m study reach, and a 3 km reach length channel slope was measured using SRTM elevation data on QGIS 3.12 to be 0.129% for the section of the Wraysbury River adjacent to Staines Moor. Median discharge is 0.360 m 3 s À1 Samples were collected from within a quadrat (basal area 0.5 m Â 0.5 m; 0.25 m À2 ; depth 25 mm) to obtain the surface grains and quagga mussel druses present on the riverbed at 20 systematically placed sample locations. These were separated by $10 m longitudinally and by $1 m laterally along a 60 m reach of channel at the study site. Within each quadrat, the b-axes of 50 randomly selected loose grains (not attached to quagga mussel druses) were measured to estimate local surface grain size percentiles and determine the size distribution of sediment available to quagga mussels. The total number of quagga mussels in the quadrat was then recorded. Twenty druses were randomly selected, and all grains > 2 mm attached to each druse (i.e., by the byssal threads of mussels) were then manually removed.
For each selected druse, the b-axis of each grain was measured and the total mass of all grains per druse weighed using portable field scales (accurate to 0.01 g). The wet mass of each quagga mussel was also measured.
The total mass of sediment attached in quagga mussel druses, and the ratio of quagga mussel biomass to mineral sediment mass was calculated. For each sample, the mean druse mass was multiplied by the observed druse density to estimate total druse mass in each sample (in g m À2 ). Materials within the 20 selected druses were separated F I G U R E 2 Location of field study site, within the context of (a) the River Thames, UK, and (b) Staines-upon-Thames [Color figure can be viewed at wileyonlinelibrary.com] into mineral sediment and quagga mussel mass. For each druse, the ratio of sediment mass to quagga mass (SMR) was calculated as: where M DS is the mass of druse sediment (in grams) and M DM is the mass of druse mussels (in grams).
The total mass of all druses (M D ) in each quadrat was then estimated as: where M M is individual mussel mass, and D M is mussel density (in m À2 ).
To investigate whether the mass of sediment attached to quagga mussel druses was dependent on the mass of quagga mussels in the druse, linear regression analysis was undertaken between the mass of mineral sediment (dependent variable) and the mass of quagga mussels (independent variable) for all of the druses from all 20 sampling locations. This was considered for (1) all druses, (2) druses containing only one quagga mussel, and (3) druses with more than one quagga mussel. To investigate whether quagga mussel density was a limiting factor for the mass of sediment attached to druses across the 60 m reach, linear regression analysis between the density of mussels in each sample (independent variable; in m À2 ) and sample SMR (dependent variable) was also undertaken.
To identify whether quagga mussels attach to all available bed sediment size fractions across the 60 m reach, regression analysis was undertaken of the D 16 , D 50 , D 84 , and sorting coefficient of the riverbed grains (independent variables) with the D 16 , D 50 , D 84 , and sorting coefficient of grains attached to quagga mussel druses (dependent variables) at each quadrat. To identify whether quagga mussels attach to all available bed sediment size fractions at the site scale, paired ttests were also undertaken comparing the difference in D 16 , D 50 , D 84 , and sorting coefficient of the riverbed grains with the grains attached to quagga mussel druses per quadrat.
All data showed a normal distribution (Shapiro-Wilk test, p > 0.05 in all cases), and linearity was assessed visually by inspection of correlation plots. Analyses were undertaken using MS Excel (Microsoft Office, 2013) and SPSS Statistics (version 27; IBM, 2020).

| Laboratory set up
Experiments were undertaken in a 10 m long, 0.6 m wide recirculating flume in the Loughborough University River Science Laboratory, UK. The flume set up consisted of roughness boards to create turbulent flow (river gravels sieved to the same size distribution as in the Wraysbury River fixed to plywood boards, raised above the flume base by 100 mm), an experimental area (460 mm Â 400 mm Â 80 mm) comprised of loose gravel and quagga mussel druses, and a sediment trap, to collect gravels mobilized from the experimental area ( Figure 3). Flume slope was set to 0.004, which was similar to the slope measured at the study site (0.0047). Experiments were not hydraulically scaled (i.e., a 1:1 scale was employed) to aid the interpretation of critical shear stress values obtained, and to increase the applicability of observed results to sediment transport modelling and to real world settings (Baynes et al., 2018).

| Experimental treatments
Three experimental treatments of different quagga mussel densities were tested; (1) a control in the absence of mussels (0 mussels m À2 ), (2) the density of mussels present in the Wraysbury River (determined by a preliminary survey; 135 mussels m À2 ), and (3) a likely invasion scenario (270 mussels m À2 ). Quagga mussel density at the time of collection (135 mussels m À2 ; August 2021) had increased by 150% since 2016 (mean 54 mussels m À2 ; Mills et al., 2017), with invasive population densities frequently exceeding 1000 mussels m À2 (Nalepa et al., 2009;Strayer et al., 1999). Therefore, 270 mussels m À2 was an appropriate density to test, as it represents a likely invasion scenario. This density was also exceeded locally at some sampling locations. Ten experimental runs were undertaken per treatment, giving a total n of 30.
All quagga mussel druses for use in the experiments were collected from the Wraysbury River in August 2021, 1 week prior to experiments (n = 750), where they were preserved in industrial methylated spirit. Dead quagga mussels were used for biosecurity reasons, but the keratinous byssal threads produced by Dreissenid mussels can preserve in the environment for long after the mussel has died (Burlakova et al., 2000), and so it was considered unlikely that preservation would significantly impact experimental results. All quagga mussels had a shell length > 25 mm, and were attached to at least two gravel clasts.
Experimental trays were filled with gravel of the same size distribution as in the Wraysbury River so that the experimental gravel surface was the same as that of the roughness boards. For treatments including quagga mussel druses, individual druses were placed manually into the experimental tray, equally spaced across the bed. The byssus-substrate components of each druse were gently worked into the top surface of the test bed while the mussel shell protruded from the surface, analogous to our observations in the field ( Figure 3).
Experimental trays were water-worked for 3 h prior to the start of each experimental run to simulate natural riverbed particle sorting processes at a depth-averaged velocity of 1.0 m s À1 and a depth of 0.25 m. The velocity was 0.1 m s À1 higher during water working than the threshold for incipient motion (in the absence of mussels, following Mills, 2019). Any sediment transported out of the test bed was fed back into the flume upstream of the test bed.

| Measurement and determination of discharge steps
Each gravel bed was subjected to series of 10 discharge steps with increasing discharge, flow velocity and bed shear stress, to examine the flow velocity and shear stress required for entrainment to occur, with the erosion of sediment measured during each step (following Mason et al., 2022). During each discharge step, flow velocity was recorded at 0.4d (where d is water depth) to estimate depth averaged velocity (U d ) using a Valeport electromagnetic flow meter. A Nortek acoustic doppler velocimeter (ADV) was used to collect velocity time series data at 25 mm above the bed. ADV data were collected for 3 min with six replicates per discharge step. ADV data were postprocessed using the velocity signal analyser tool (Jesson et al., 2015) and underwent phase-space thresholding with spikes replaced by linear interpolation (Biron et al., 2004;Goring & Nikora, 2002). Data with correlation below 70% or signal-to-noise ratio below 15 were removed (following McLelland & Nicholas, 2000;Nortek, 2009).
Velocity time series data were used to quantify near bed velocity (U b ) and bed shear stress (τ b ) using the turbulent kinetic energy (TKE) approach.
TKE (E) was calculated according to: where u 02 , v 02 and w 02 are the second order moment statistics for velocity in the stream-wise, cross-stream and vertical directions, respectively. Bed shear stress τ b (in N m À2 ) was then determined as: where C 1 is a constant (the value used was 0.19; Soulsby, 1983;Stapleton & Huntley, 1995), and ρ w is water density (1000 kg m À3 ).
The values U d , U b and τ b were tested for differences between discharge steps using analysis of variance (ANOVA) and post hoc Tukey honest significant difference (HSD) tests. All discharge step treatments were different, with the exception of near bed shear stress and near bed flow velocity between step 3 and step 4. (p < 0.05). Discharge steps therefore ranged from a flow velocity and shear stress below the motion of any particles up to the maximum velocity and shear stress achievable in the flume in the described set up, and well above entrainment thresholds (range in velocity was 0.4 to 1.3 m s À1 at 60% depth, increasing in 0.1 m s À1 increments; range in shear stress 0.52 to 5.55 N m À2 ).

| Experimental procedure
The flume was filled to a depth of 0.25 m above the fixed boards and experimental area, and the pump speed was increased to the first discharge step. Each discharge step lasted for 180 s. At the end of the discharge step, any sediment caught in the bedload sediment trap was removed using a 1 mm mesh net. The pump speed and weir height were then adjusted to increase the flow to the next discharge step, and the process was repeated through 10 discharge steps.
The sediment removed from the sediment trap at the end of each discharge step was sieved into half-phi grain size fractions. This allowed for the total mass of sediment and the mass of sediment of each grain size at each discharge step to be analysed. The total mass of any quagga mussel druses recovered from the trap at each discharge step was also recorded.

| Data analyses
We defined incipient motion of sediment as 1% of surface material comparisons; here and throughout) were undertaken to examine differences in τ c of the bed between experimental runs with differing quagga mussel densities.
The Shields parameter (τ *c ; Shields, 1936) was calculated for each treatment because it provides a useful means to compare critical entrainment thresholds across grain sizes and is frequently used in modelling bedload sediment transport (e.g., Meyer-Peter & Müller, 1948).
where τ c is critical shear stress, ρ s is sediment density (2650 kg m À3 ), g is gravitational acceleration (9.81 m s À2 ), and D 50 is sediment grain size (D 50 = 12.4 mm).
This was plotted against grain Reynolds number (Re * ) to compare the entrainment thresholds of each treatment to the curve developed by Shields (1936), and to other stabilizing taxa reviewed by Mason and Sanders (2021).
To assess if quagga mussels affect sediment transport after the threshold for incipient motion has been reached (bedload flux), the cumulative mass of sediment mobilized across each experimental run, and the mass of sediment mobilized at each discharge step, were compared between mussel density treatments using Kruskal-Wallis (H) tests.
Whether quagga mussels affect the mobility of different grain sizes was also analysed. A flood competence function (cf., Montgomery & Buffington, 1997) was used to understand how quagga mussels affected the τ c of each grain size of sediment; Kruskal-Wallis (H) tests were used to compare the shear stress at which the movement of each grain size was first detected. Finally, to analyse whether quagga mussels affect the grain size distribution of mobilized gravels, the D 16 , D 50 , and D 84 of the grains mobilized at each discharge step were calculated. Kruskal-Wallis (H) tests were used to analyse differences in D 16 , D 50 , and D 84 between treatments with different quagga mussel densities at each discharge step, with Spearman's Rank correlation tests used to analyse if the D 16 , D 50 , and D 84 changed with increasing discharge steps.

| Bedload sediment mobility estimates
To increase the scale of our study and its application to river geomorphological and invasive species management, we modelled the number of days in which sediment transport is likely to occur in the presence of 0, 135, and 270 quagga mussels m À2 (Q4). Daily discharge data from the Wraysbury River were used to model the reach averaged boundary shear stress (τ 0 ) and thus partitioned grain stress (τ g ) time series. We then calculated the proportion of days that the critical shear stress (and thus critical grain stress) observed in flume experiments exceeded those modelled for the Wraysbury River.
It is important to recognize limitations to this method; the slope and depth of the flume channel differ from that of the Wraysbury River, the calculation of τ g is different for the flume channel (TKE from ADV measurements) and the river channel (Manning's equation), and τ is recognized to be a poor descriptor of bedload sediment transport because it does not consider spatial and temporal variability in bedload transport which is typically high . Thus, the goal of this modelling is not to provide explicit sediment flux predictions, nor to reconstruct bedload sediment transport over the period of quagga mussel invasion. Rather, this modelling exercise aimed to assess whether changes to critical shear stresses enacted by the presence of quagga mussels are of sufficient magnitude to have important geomorphic implications to the mobility of sediment in a river.
Reach averaged boundary shear stress was then estimated as: where τ 0 is reach averaged boundary shear stress (in N m À2 ), ρ is water density (1000 kg m À3 ), g is gravitational acceleration (9.81 m s À2 ), and S is river slope.
However, only a portion of the boundary shear stress is available for bedload transport. Therefore, τ 0 was partitioned to grain stress (τ g ) according to Wilcock et al.'s (2009) application of the Manning-Strickler formula: where n is Manning's n (n = 0.045, commensurate with gravel bed rivers including the Wraysbury River), and n D is the Strickler relation (n D = 0.020 for measured study site D 65 of 14.3 mm).
Applying n = 0.045, commensurate with gravel bed rivers including the Wraysbury River, and using D 65 for D (as in Wilcock et al., 2009;D 65 for study site = 14.3 mm), the relationship between τ 0 and τ g at the Wraysbury River can be expressed as: This relationship between τ 0 and τ g was applied to the estimated τ 0 from daily discharge data to obtain a daily grain stress series, and to the critical shear stress results observed in flume experiments to estimate critical grain stress. We then calculated the number of days when the critical grain stresses recorded in the flume in the presence of 0, 135, and 270 mussels m À2 were exceeded in the field between 2014 and 2019. A comparison of these values allowed for an estimation of the number of days per year where the incipient motion of gravel was inhibited by the presence of quagga mussels.
There was a significant but weak positive relationship between the mass of quagga mussels and the mass of mineral grains in a druse ( Figure 4a). This relationship was slightly stronger when only druses with more than one quagga mussel were considered (Figure 4b), and absent when druses with only one mussel were considered (r 2 = 0.002, p = 0.429; Figure 4c). The mass of sediment attached to druses in each sample appeared to be dependent on local quagga mussel abundance; as quagga mussel density increased, the ratio of mineral sediment to mussel mass in a druse (SMR) reduced ( Figure 4d).
At the reach scale, the D 16 , D 50 , D 84 , and sorting of the quagga mussel druses was dependent on the D 16 , D 50 , D 84 , and sorting of the F I G U R E 4 Individual druse sediment and quagga characteristics. Relationships between quagga mussel mass and grain mass in druses when (a) all druses, (b) only druses with one mussel, (c) only druses with more than one mussel are considered, and (d) between quagga mussel density and the sediment to mussel mass ratio. Asterisks indicate significant correlations (*p < 0.05; **p < 0.01; ***p < 0.001) bed grain material across all quadrats. (Figure 5a,b). At the quadrat scale, the grains attached to quagga mussel druses were significantly finer (D 16 = À2.6 mm; D 50 = À4.0 mm; D 84 = À3.4 mm) and more poorly sorted (+0.14 mm) than the sampled surface grains across all samples (Figure 5c).

| Flume experiments
Quagga mussels significantly increased the critical shear stress for incipient motion of the riverbed (τ c ; H (2) = 17.326, p < 0.001). Critical shear stress for the control treatment gravels was 2.94 N m À2 , which increased by 20% to 3.52 N m À2 , and by 40% to 4.12 N m À2 in the presence of the medium and high-density quagga treatments, respectively. This corresponded to an increase in Shields criterion from a value of 0.015 in the control treatment to 0.018 and 0.021 in the medium and high-density treatments, respectively ( Figure 6). The difference in τ c between the zero-density treatment and the medium density treatment was not significant (p = 0.382), but τ c was significantly higher for the high-density treatment than the medium (p = 0.029) and zero (p < 0.001) density treatments (Figure 7a). This was also observed considering specific grain sizes, with both medium and high-density quagga mussel treatments increasing mean observed τ c across all grain fractions (Figure 7b).
The cumulative mass of sediment mobilized was significantly lower in the high-density treatment than when quagga mussels were absent, for every discharge step (step 6, p = 0.004; steps 7-10, p < 0.001; Figure 8a), and significantly less than the medium density treatment in the higher discharge steps (step 8, p = 0.029, steps 9 and 10, p < 0.01). The cumulative mass of sediment mobilized in the medium quagga mussel density treatment was not different to when quagga mussels were absent, with the exception of discharge step 5 (p = 0.003; steps 6-10, p > 0.05).
A similar pattern was observed when the mass of sediment mobilized at each discharge step was considered (Figure 8b). A lower mass of sediment was mobilized in the high discharge steps in the high-F I G U R E 5 Differences and relationships between riverbed sediment and quagga mussel druse sediment. (a) The relationships between the D 16 , D 50 , and D 84 of riverbed and druse sediment; (b) the relationship between the sorting of riverbed and druse sediment; and (c) differences between the D 16 , D 50 , D 84 , and sorting of riverbed and druse sediment. Asterisks indicate significant correlations, and significant pairwise comparisons (*p < 0.05; **p < 0.01; ***p < 0.001) F I G U R E 6 Effects of biostabilization on Shields criteria, considering quagga mussel (Dreissena bugensis) druses in this study, and previous studies reporting changes to entrainment thresholds of river gravels by net spinning caddis larvae (Albertson et al., 2014(Albertson et al., , 2019Johnson et al., 2009), following Mason and Sanders (2021). Arrows show approximate change in Shields parameter and Reynolds number which occurs due to the biostabilization of sediment [Color figure can be viewed at wileyonlinelibrary.com] density treatment than the medium density (steps 9 and 10, p < 0.01) and control treatment (steps 7-10, p < 0.01). The mass of sediment mobilized at each discharge step in the medium quagga mussel density treatment was not different to when quagga mussels were absent, with the exception of step 5 (p = 0.003; steps 6-10, p > 0.05).
The grain size distribution of mobilized material was generally consistent throughout all runs (Figure 8c). Mobilized grain size was predominantly independent of critical shear stress, and the D 16 , D 50 , and D 84 of the mobilized sediment were generally consistent between mussel densities and at all discharge steps.

| Bedload sediment mobility estimates
In the absence of quagga mussels, incipient motion was estimated to occur on 29.4% of days of the 5-year time series (discharge ≥ Q 30 ).
The number of days on which incipient motion was modelled to occur reduced substantially in the presence of quagga mussels. Incipient motion was modelled to occur on just 6.5% (≥ Q 7 ; medium density) and 1.5% (≥ Q 2 ; high density) of days ( Figure 9). During the 5-year time series, there were only 10 separate (> 48 h separation) flow events modelled to exceed the shear stresses required for incipient motion in the presence of the high-density quagga mussel treatment.

| DISCUSSION
Quagga mussels are among the most invasive species in the world (Roy et al., 2014), and this study demonstrates that invasion by quagga mussels may lead to substantial change to sediment mobility.
Quagga mussel druses attached to a large mass of riverbed particles, and to all sizes of riverbed particles available. Quagga mussel attachment resulted in increased shear stresses required to entrain bed gravels, and is predicted to reduce movement of riverbed grains during high flow events. Applying these shear stress values to high resolution flow data from the Wraysbury River indicated that quagga mussels may reduce the number of days per year where sediment transport occurs, with medium and high densities modelled to increase the required discharge for a geomorphically active flow event from Q 29 to Q 7 and Q 2 , respectively. These results reinforce the current evidence that invasive species are important drivers of geomorphic change (Crooks, 2002;Emery-Butcher et al., 2020;Fei et al., 2014;Harvey et al., 2011;Mason & Sanders, 2021), and extend the scope of observed zoogeomorphic riverbed stabilization, knowledge of which is predominantly limited to net building caddisfly larvae (Mason & Sanders, 2021). Further, the results presented here extend previous research by presenting the first multi-year modelling of zoogeomorphic changes to sediment mobility due to stabilizing taxa.
4.1 | Quagga mussels are not selective of grain size Quagga mussels attached together riverbed sediment (mean 1114 g m À2 inclusive of mussels, 591 g m À2 of exclusively mineral sediment). This is a greater mass than has previously been recorded for other taxa that attach particles with secretions (e.g., case building caddisfly larvae; mean = 38 g m À2 , maximum = 139 g m À2 ; Mason et al., 2019; $ 300-470 g m À2 Statzner et al., 2005) due to the larger size of particles used by quagga mussels. Whilst many case-building caddisfly species are selective in the grain sizes used to construct cases (Mason et al., 2019), this study suggests that quagga mussels will attach to all available riverbed sediment. Reach scale analysis of quagga mussel druses demonstrated a significant relationship between the sediment size distributions of grains attached to druses and bed grains, suggesting that druse attachment was driven by grain availability, rather than active selection. Whilst it is not possible to quantitatively extrapolate beyond the sediment grain sizes measured here (D 50 of 9 mm to 19 mm at the site scale), quagga mussels attaching to all available grain sizes suggests that similar effects to those observed in this study are also likely to be observed in rivers with coarser or finer grained gravel beds than examined here.
At the site scale, the grain sizes attached to druses were significantly finer than the surface bed material, but the magnitude of this difference was small (D 16 = À2.6 mm; D 50 = À4.0 mm; D 84 = À3.4 mm). The reason for this difference in grain size may reflect the sampling strategy employed. Quagga mussels did not omit to attach to coarse sediment, rather they also attached to subsurface fine sediment particles, which were not accounted for in sampling surface material.
F I G U R E 7 Critical shear stress required (a) for incipient motion of the riverbed in the presence of 0, 135, and 270 quagga mussels m À2 , and (b) for incipient motion of each grain size fraction (finer than) in the presence of 0, 135, and 270 quagga mussels m À2 . Asterisks indicate significant pairwise differences between mussel density treatments (*p < 0.05; **p < 0.01; ***p < 0.001) Another control on the mass and size of sediment that quagga mussels attached to was the density of quagga mussels in the sample.
As quagga mussel density increased, the mass of the individual druses decreased. This could suggest a density dependent factor; that all available coarse sediment had been 'used'. However, high numbers of dead quagga mussel shells (which have a low mass density) were found present on the substrate during our field study, and were attached to many of these druses. Bivalve shells can persist for hundreds to thousands of years in fluvial environments (e.g., Himson et al., 2020;Huntley et al., 2021;Wolverton et al., 2010), and so it is reasonable to assume that across decades, new populations of mussels may be forced to establish on beds of dead mussels from prior generations, resulting in a lower mass of abiotic grains attached to the sampled surface druses. Multiple layers of quagga mussels would likely give increased protection to bed gravels and further reduce sediment transport, with surface mineral grains buried under multiple layers of quagga mussels. This is particularly likely given that byssal threads created by mussels can be preserved in the environment long after the mussels have died (Burlakova et al., 2000), suggesting that dead mussels attached to the grains will still have a geomorphic effect via attachment.

| Quagga mussels increase critical entrainment shear stress
The attachment of grains by quagga mussels significantly increased the threshold for the critical entrainment of gravel in flume experiments. The results of this study are consistent with those of Mills (2019), which observed that a medium quagga mussel density of 125 mussels m À2 did not significantly reduce sediment transport during a high flow velocity that exceeded the critical entrainment threshold, but that a high density of 250 mussels m À2 did. Whilst Mills (2000 m À2 ) of silk spinning caddisflies (0.14-0.92). This indicates that the densities of quagga mussels studied here (135 and 270 m À2 ), which are low compared to densities of many invaded sites (> 1000 individuals m À2 ; Nalepa et al., 2009;Strayer et al., 1999), may stabilize gravels to a similar extent to that of natural densities of caddisflies.
Established quagga mussel populations may therefore stabilize gravels at least to the degree that has been observed for net-building caddisflies, which have the greatest known stabilizing capacity of any freshwater invertebrates (Mason & Sanders, 2021).
The biostabilization observed in these flume experiments examined quagga mussel density alone. However, total biostabilization may be enhanced by the presence of other animals. The abundance of net building caddisfly was positively associated with quagga mussel density in the Wraysbury River (Hydropsyche spp.; Mills et al., 2017). This effect may also be non-linear; Albertson et al. (2014) found that at high densities, the stabilization of river gravels by polycultures of multiple caddisfly species was greater than the expected additive effects F I G U R E 8 Mass of sediment mobilized in the presence of 0, 135, and 270 quagga mussels m À2 , considering (a) cumulative mass, and (b) mass mobilized in each discharge step, and (c) mean sediment grain size (AE1 standard error) of sediment mobilized during each run in the presence of 0, 135, and 270 quagga mussels m À2 . Asterisks in (b) indicate significant pairwise differences between mussel density treatments (*p < 0.05; **p < 0.01; ***p < 0.001) of the independently tested species. Understanding the total community level biostabilization of sediments is therefore an important avenue for future research to understand the full extent of quagga mussel biostabilization.
4.3 | Quagga mussels potentially reduce annual bedload sediment flux Applying the critical grain stress results observed in the flume experiments to the estimated grain stress series in the Wraysbury River indicated that the presence of quagga mussels reduced the number of days on which incipient motion was modelled to occur, from 29.4% of days in the absence of quagga mussels to just 6.5% and 1.5% of days in the presence of 135 and 270 mussels m À2 , respectively. To examine what affect this may have on bedload sediment flux, we use the Wong and Parker (2006) version of the Meyer-Peter and Müller (1948) bedload transport equation with discharge time series for water years 2014-2019 to estimate cumulative bedload yield (cf., Albertson et al., 2014): where q* b is volumetric unit bedload transport (in kg s À1 ), W is river width (6.4 m), and r = (ρ s -ρ w )/ρ w ., τ *c was applied from the flume experiments, and τ c was calculated for the study site as: where D = 12.4 mm and S = 0.00129.
We present this as a discussion point rather than as a result, because of known uncertainties in the application of Meyer-Peter and Müller (1948) equations (e.g., see Gomez & Church, 1989;Hinton et al., 2018;Wilcock et al., 2009), and because this modelling represents a conceptual extrapolation from laboratory flume results to estimate indicative proportional changes under current and future invasion scenarios (sensu Albertson et al., 2014).
For the full time series, the medium and high densities of quagga mussels reduced total estimated bedload yield by 57.0% and 74.0%, respectively, compared to when quagga mussels were absent. This varied by year (Figure 10), with the majority of sediment moved in the presence of quagga mussels during flood events in 2014-2015 and 2015-2016. In the 3 years where flooding did not occur on this scale, the medium and high densities of quagga mussels reduced total estimated bedload yield by > 81% and > 95% of bedload sediment, respectively.
This modelling exercise represents a hypothesis of the potential effects to bedload transport that quagga mussel druses may have in the field. Field validation of this modelling using methods such as bedload sampling during high flow events (cf., Laronne et al., 2003;Rickenmann & McArdell, 2007), tracer particles to evaluate sediment transport during high flow events (cf., Hassan & Roy, 2016), field flumes to accurately observe the processes in situ during a high flow event (cf., Batalla et al., 2010), and field surveying before and after a high flow event would aid in the validation of this conceptual modelling. Nevertheless, this exercise indicates the potential for alteration to fluvial sediment budgets following quagga mussel invasion.
4.4 | Quagga mussels and their effects are geographically widespread, and have ecosystem-wide importance The stabilization of gravel beds during high flow events may have important ecological and biogeochemical implications beyond changes to channel sediment dynamics and morphology. The biostabilization of sediments creates stable habitat patches which may act as refugia for invertebrates during high flow events (Cardinale et al., 2004;Lancaster, 2000). However, reduced gravel movement may also result in increased rates of fine sediment clogging, as high flows may be less able to rework surface material and reduce colmation (Mathers et al., 2021;McKenzie et al., 2022). This can alter stream bed processes through many mechanisms (reviewed in Wharton et al., 2017), with potentially deleterious impacts for macrophytes (Jones, Murphy et al., 2012), diatoms (Jones et al., 2014), macroinvertebrates , and fish (Kemp et al., 2011). High levels of fine sediment were observed in the reach in areas of high quagga mussel density, which may impact biology in the Wraysbury River.
Further, the stabilization of gravel beds can initiate a positive feedback cycle; reduced bedload transport can increase downstream bed armouring (Piqué et al., 2017), leading to increased stability against high flow events, and a subsequently reduced bedload sediment flux downstream of the initially stabilized area (Brenna et al., 2021). Therefore, by starving downstream reaches of sediment, quagga mussels have the capacity to change the bedload flux beyond invaded reaches.
In this respect, future study will be important. The stabilizing effects of quagga mussel observed in this study are likely to become more pronounced in space and in time. Following the detection of quagga mussels in the Wraysbury River in 2014 , past studies of quagga mussels on the Wraysbury River measured a F I G U R E 9 Modelled bed grain stress in the Wraysbury River. (a) Modelled reach averaged grain stress (τ g ) at the Wraysbury River. Horizontal lines denote the riverbed τ g from flume experiments in the presence of 0, 135, and 270 quagga mussels m À2 . Modelled shear stress exceeding the τ g lines indicates a geomorphically active flow event mean mussel density of 54 mussels m À2 from 2015 to 2016 (Mills et al., 2017), and this study observed mean densities of 122 mussels m À2 (2021). Whilst the densities of quagga mussels observed in this study indicate a doubling of the population in 5 years, such densities are still low compared to other invaded environments, where > 1000 individuals m À2 have been frequently observed (Nalepa et al., 2009;Strayer et al., 1999). These population densities may be achieved rapidly in invaded catchments. For example, Dreissena polymorpha densities increased from 173 to 8816 mussels m À2 in Lake Michigan over just a 10-year period (Nalepa et al., 2009). Of 591 potentially invasive taxa, quagga mussels were identified as the invasive species posing the highest risk to invade, establish, and have negative effects throughout the UK (Roy et al., 2014). This risk may be similar globally; quagga mussels have been widely reported to be replacing the closely related and more widespread zebra mussel (D. polymorpha) throughout North America (Matthews et al., 2015;Ram et al., 2012). Continued temporal monitoring of sites such as the Wraysbury River is needed to understand the potential magnitude of invasion, and thus geomorphic change, that quagga mussels could present.
As well as increasing in density, Dreissenid mussels have also been found to rapidly expand spatially throughout freshwater systems elsewhere (mean 120 km a À1 northwest Europe; Matthews et al., 2014), meaning that the stabilization of riverbed substrate following quagga mussel invasion may become a widespread effect.
Indeed, these effects are to be potentially observed on a global scale; as of 2021, quagga mussels have been recorded as present in at least 57 locations throughout the southeast of England (NBN Atlas, 2021a; Figure 11b), in at least 10 countries, and at least 16 states in the United States (CABI, 2021a; Figure 11a). Further, the closely related and more widely distributed zebra mussel also uses byssal threads to anchor itself to aquatic substrates, and thus may have similar stabilizing effects to quagga mussels (Mason & Sanders, 2021). Zebra mussels are present in at least 1405 locations throughout the UK (NBN Atlas, 2021b; Figure 11c), in at least 36 countries, and at least 28 states in the United States (CABI, 2021b; Figure 11a), and so their potential geomorphic effects in rivers are worthy of study.
4.5 | Invasive species are particularly potent geomorphic agents The study of biogeomorphic processes associated with invasive taxa is important. Whilst this study has suggested that quagga mussels may present a similar biostabilizing capacity as populations of native F I G U R E 1 0 Modelled bedload yield in the study section of the Wraysbury River. Panels (a) to (e) show estimated cumulative bedload transport for each hydrological year, modelled using Wong and Parker's (2006) amendments to Meyer-Peter and Müller's (1948) bedload transport equations. Values in the top right corner of each panel indicate the reduction in bedload yield associated with the presence of the medium and high density of quagga mussels. 'Water day' refers to days relative to the start of the hydrological year (1 October) [Color figure can be viewed at wileyonlinelibrary.com] caddisfly larvae, invasive species have the potential to cause stronger geomorphic disturbances than native species. The activities of native taxa are likely reflected in long-term landscape evolution, due to the temporal persistence of these processes (e.g., Fremier et al., 2018), whereby landscapes have reached an equilibrium with biogeomorphic processes in place. For example, a river with a slope and sediment grain size already adjusted to the behaviours of native species may return a relatively low biotic component of the sediment budget, as their effects have already caused system wide adjustment. However, in landscapes that have formed independently to the presence of the invasive species, the abrupt change brought about by an invasive species may force rivers into a status of non-equilibrium. Thus, invasive species are able to enact rapid threshold changes within invaded systems due to the disparity of the landform, and the abruptly altered process in the presence of the invader.

| CONCLUSION
The invasion of quagga mussels will likely present abrupt changes to geomorphic processes in invaded catchments. In the Wraysbury River, quagga mussels attached to a high mass of surface grains, which comprised all size fractions of available sediment. Population densities representing a likely invasion scenario, which were also found locally, significantly increased the critical shear stress required for the incipient motion of riverbed gravels in flume experiments by 40%. Further, these were estimated to increase the discharge requirement for a geomorphically active flood event from Q 30 to Q 7 and Q 2 . These combined results demonstrate that quagga mussels have the capacity to alter fluvial sediment transport processes, which may have important effects for channel morphology, ecological functioning, and human infrastructure in invaded systems. Whilst this study has focused on one study site, these effects and processes are likely to become more widespread, due to the increasing distribution and highly invasive nature of the quagga mussel.
Aquatic biosecurity from invasive species therefore is a key geomorphological, as well as ecological, management consideration.
Whilst some sediment transport and stream evolution models (e.g., Castro & Thorne, 2019) are starting to conceptually incorporate biotic energy, biologyand in particular invasive speciesis still largely disregarded . As such, future geomorphic modelling and management must seek to include biotic energy as a key component to fully understand the rates and processes of landscape development.
F I G U R E 1 1 Distribution of quagga mussel and the closely related zebra mussel (Dreissena polymorpha). (a) Global distribution of Dreissenid mussels. Data display the presence of mussels recorded by country, except for the United States and Canada, which are mapped by state and province, respectively, due to the availability of data. Data from CABI (2021a, 2021b). The native ranges of Dreissena bugensis and D. polymorpha are restricted to a small number of estuaries, deltas, and lowland river channels in the northern Black Sea region (described in Son, 2007;Zhulidov et al., 2010