A numerical model for simulation of freshwater mussel dynamics was used to investigate the effects of substrate and hydrodynamic conditions on the formation of mussel beds in a 10-km reach of the Upper Mississippi River (UMR). Suitable habitats for mussel survival were identified by creating a dimensionless parameter (shear stress ratio) combining shear force and substrate type. This parameter is a measure of substrate stability that could be used in many different applications. Dispersal of juvenile mussels with flow as they detach from their fish hosts was simulated by a particle-tracking mechanism that identified suitable areas for colonization with the potential to evolve into mussel beds. Simulated areas of mussel accumulation coincided with reported locations of mussel beds, and simulated densities were in the range of abundant mussel beds in other reaches of the UMR. These results, although more qualitative than quantitative, provide insight into factors influencing the formation of mussel beds in a large river.
Mussels of the Family Unionidae are worldwide in distribution, but their maximum diversity and abundance occurs in North America where 297 species have been documented (Turgeon et al. 1998). Unionids probably provide critical ecosystem functions such as particle processing, nutrient release, and sediment mixing, but these functions have not been thoroughly addressed (Strayer et al. 2004, Vaughn et al. 2004). The combined effects of overharvesting, habitat alteration, and invasive species have placed unionids among the most endangered faunal groups (Strayer et al. 2004).
Unionids are long-lived animals (30–100 y; Bauer and Wachtler 2001) with a complex life cycle. They have specialized larvae (glochidia) that are obligate parasites on fishes. Glochidia remain attached to a fish host for a period of weeks to months (species and temperature dependent) while they undergo organogenesis. In this way, they may be dispersed into different habitats (Vaughn and Taylor 2000). Juvenile mussels detach from the fish host to begin their lives as free-living organisms. Successful reproduction requires the availability of an appropriate host fish at the appropriate time and that juveniles find favorable habitat. The survival of mussels is dependent ultimately on the interaction of several biotic and abiotic factors operating at different spatial and temporal scales (Haag and Warren 1998, Vaughn and Taylor 2000, Strayer et al. 2004).
Unionids are patchily distributed and exist in aggregated multispecies groups called mussel beds (Strayer et al. 1994). Mussel beds may be the product of differential mortality in which juveniles settle evenly on the river bottom but perish in unsuitable habitats (Strayer 1999). Conversely, areas devoid of juveniles could simply mean that juveniles have not yet arrived there (Vaughn and Taylor 2000). Attempts to discriminate micro- and macrohabitat differences between areas with and without unionids have been largely unsuccessful (Holland-Bartels 1990, Haag and Warren 1998, Strayer 1999), but several hydraulic variables recently have shown some success in predicting the distribution of unionids (Layzer and Madison 1995, Hardison and Layzer 2001, Howard and Cuffey 2003). When taken together, these data suggest that conditions that reduce the shear forces acting on river bottoms may be important in creating refuges from flow that enable unionids to persist under high-flow conditions.
Unionids live partially or completely buried in the sediments of rivers, so substrate and hydrodynamic conditions probably have an effect on community structure. Mussels require appropriate substrate to anchor and burrow, and nearby bed currents, mean velocities, and substrate stability also may affect their distribution (Holland-Bartels 1990, Layzer and Madison 1995, Strayer 1999, Hardison and Layzer 2001, McRae et al. 2004). Interactions among this suite of physical factors make their individual effects difficult to distinguish (Rempel et al. 2000, Arbuckle and Downing 2002, Jowett 2003). In addition, many traditional variables (i.e., water depth and flow) are flow conditional and measures made at one discharge are of limited value in predicting habitat suitability at other flows (Layzer and Madison 1995, Strayer 1999).
Our paper reports on the application of a numerical model to analyze the effects of substrate and hydrodynamic conditions on habitat suitability and the formation of mussel beds. A dimensionless parameter for the assessment of substrate stability was created (shear stress ratio) and tested against published data on habitat preferences for mussels. This parameter was combined with reported data on preferences of unionids related to traditional hydrodynamic variables (water depth, flow velocity, and substrate type) to estimate habitat suitability under various flow regimes. The life cycle of freshwater mussels was simulated with a population dynamics model representing the various life stages, their dependence on appropriate habitat conditions for survival, the interaction between glochidia and host fish, and the dispersal dynamics of juvenile mussels. The model was applied to a 10-km reach of the Upper Mississippi River (UMR) and the simulated areas of mussel accumulation were compared with historic data on the actual location of mussel beds.
Substrate stability assessment
A nondimensional parameter, the shear stress ratio (RSS), was created as a measure of substrate stability. The RSS is defined as:
The value of τc is controlled by shear stress, flow velocity, riverbed materials, and armoring (Lorang and Hauer 2003); unfortunately, no simple expression exists to predict τc over a broad range of hydrodynamic conditions and sediment types. τc was defined in our work according to the Shields diagram (Shields 1936); nevertheless, other equations can be used (e.g., van Rijn 1993, Wilcock 1993, Lorang and Hauer 2003).
High flows are the limiting conditions for substrate stability and, therefore, τ0 and RSS must be computed at appropriate flow rates. Strayer (1999) noted that measures at low flows may not be obviously related to the location of mussel beds and suggested that mussel beds may generally be found in areas that remain stable during flows with recurrence intervals between 3 and 30 y. Given the difficulties of gathering field measurements at high flows, a numerical model becomes an extremely useful tool for estimating τ0 at high flows.
We hypothesized that higher mussel densities would be found where RSS < 1 under most flow regimes. However, a maximum tolerance of RSS = 1.25 was set because some species may withstand significant sediment motion (Di Maio and Corkum 1995). This concept was tested against data by Layzer and Madison (1995), and our analyses show that reported water depth and substrate preferences for unionids can be explained largely in terms of substrate stability (Appendix).
Habitat suitability assessment
A habitat suitability index model was developed for identification of potential mussel beds based on substrate and hydrodynamic variables (Fig. 1). Representative habitat preferences for unionids were based on literature-derived values for water depth (Miller et al. 1987, McMahon and Bogan 2001 and references therein), flow velocity (Miller et al. 1987, Holland-Bartels 1990, Johnson et al. 2001), and substrate composition (Stern 1983, Miller et al. 1987, Strayer and Ralley 1993, Strayer 1999). The RSS concept was added for substrate stability assessment.
Habitat suitability was evaluated for various flow rates covering the expected yearly hydrograph to identify those areas in the domain offering appropriate habitat for mussels throughout the year. At each cell, the overall habitat suitability score (HSS) was computed as the minimum suitability index value obtained throughout the year. The domain was divided into a 2-dimensional (2D) lattice for identification of suitable patches for mussel development.
Simulation of the life cycle of freshwater mussels
The mussel dynamics model developed by Morales-Chaves (2004) was used to simulate the life cycle of freshwater mussels. The life span of mussels was divided into 4 life stages and, at each stage, survival was determined by background mortality and habitat suitability. Background mortality probabilities were 0.9999, 0.60, 0.50, and 0.05 for the egg/larvae, parasitic, juvenile, and adult stages, respectively (Stein 1973, Bauer 2001, Jansen et al. 2001). Habitat suitability was based on the habitat suitability model described above. The cycle starts as the glochidia are expelled from the female. Glochidia that encounter a host fish move into the parasitic stage and re-enter the flow 25 d later as juveniles. Glochidia that do not encounter a host fish remain suspended in the flow for a maximum of 5 d, after which they die (value based on Jansen et al. 2001). Juveniles settling in suitable habitats survive and continue their development. Adulthood is reached after 5 y, and reproduction begins. Adult mussels live for 40 y and, at the end of their life span, the population continues with successive generations of mussels produced in the study area. The duration of each life stage is species- and temperature-dependent and can be modified by the user; the values presented here are generalizations based on the literature (Jansen et al. 2001, McMahon and Bogan 2001 and references therein, Haag and Staton 2003).
A key process simulated in the model is the hydrodynamic transport of larvae and juveniles as they detach from the host fish. Individual mussels were modelled as sediment particles suspended in the flow, and the settling velocity of each juvenile was computed as a function of its size and body density following principles from sediment transport (van Rijn 1993, Chang 1998). Using 0.02 cm as the mussel size (Wachtler et al. 2001) and a density of 1.01 g/cm3 (based on sinking velocities for marine larvae; Wildish and Kristmanson 1997), a settling velocity of ˜0.03 cm/s was estimated. The spatial distribution of settling juveniles was estimated using a particle-tracking mechanism used to locate a given mussel at all times.
A detailed description of the mussel dynamics model is beyond the scope of our paper. Readers interested in tracing the model formulation step by step are directed to Morales-Chaves (2004), where further details, including model parameters and values, can be found.
A 10-km river reach of Navigation Pool 16 in the UMR was chosen as a case study. Comprehensive data on the spatial distribution of substrate types in the study area were lacking, so a substrate coverage map (Fig. 2A) was created based on data from an upstream reach. Data on substrate composition in different habitat types in Pool 13 (˜75 km upstream of the study area) were available from the Long Term Resource Monitoring Program (USGS 2003); we assumed that equivalent habitat types in Pool 16 would have similar substrate composition. The main navigation channel was classified as sand according to observations by USACE (1981) in the study area. For each substrate type, a mean particle size was assigned based on typical values for the particles and in agreement with data in Holland-Bartels (1990) for an upstream reach of the UMR.
A 3-dimensional (3D) hydrodynamic model of the study area was created by the Iowa Institute of Hydraulics Research (IIHR–Hydroscience and Engineering) of the University of Iowa. The model formulation follows Lai et al. (2003). Bathymetric (cross-sectional) and velocity data used in the development and verification of the model are described in Young et al. (2005). Near-bed velocity distributions and shear stress were extracted directly from the simulation results, and a 2D depth-averaged velocity distribution was computed to characterize mean flow patterns. Six steady state scenarios were simulated (Fig. 2B) covering a range from 2 to 99% of the daily observed flows at Lock and Dam 16 and a maximum return period of 4.2 y. The 18-y flow record was used to compute the mean yearly hydrograph (Fig. 2C) used in the mussel dynamics simulations. Hydrodynamic conditions for each day of the year were interpolated linearly from the simulated flow scenarios.
Data on the spatial distribution of fishes were obtained from a resource inventory of the UMR (USACE 1984) that roughly identified spawning areas, areas used in commercial fisheries, and sport fishing areas throughout Pool 16, but did not detail habitat use by different species or temporal variation throughout the year. The locations of historic mussel beds were obtained from studies by the US Army Corps of Engineers (USACE 1981, 1984); unfortunately, these studies were largely qualitative and, thus, densities of unionids in this reach were unavailable.
A coarse grid with cell sizes varying from 10 ×10 m to 10 × 25 m was used for the mussel dynamics simulation because of the large spatial area considered in the model (˜15 km2). At this scale, it should be possible to observe general patterns of mussel-bed distribution, but not patchiness within a given bed. The spatial distribution of suitable habitats was computed for the 6 flow scenarios indicated in Fig. 2B and used as input data to the mussel dynamics model along with the hydrodynamic files describing the mean yearly hydrograph and the spatial distribution of host fish. In this simulation, fish distribution was assumed to be static because data on temporal variation in fish distributions were unavailable. The initial mussel population consisted of 10 glochidia per cell, uniformly distributed throughout the domain. No predefined information about existing mussel beds was given; instead, mussel beds were allowed to emerge according to the functional processes governing population dynamics. The model was run for 70 simulated years to observe areas of mussel accumulation and estimate mussel density.
Identification of suitable habitats
Shear forces acting on the river bottom increased as a function of increasing flow rates such that, under high flows, the onset of sediment motion was exceeded in most of the domain (Fig. 3A–C). Nevertheless, ˜38% of the stable patches observed at the average flow (2039 m3/s; Fig. 3B) remained stable even at the mean peak annual flow (3965 m3/s; Fig. 3C). These patches are areas of flow separation produced by the geometry of river banks and islands that are less affected by the flow increment than other areas. They are analogous to the flow refuges discussed by Strayer (1999) and, other conditions being appropriate, are likely to be suitable areas for mussels to inhabit.
The spatial distribution of suitable habitats varied substantially with flow rate (Fig. 4A–C). In general, as flow increased, the number of suitable areas for unionids decreased. Under low flows (e.g., 566 m3/s, the flow with 99% probability of exceedance), most of the substrate remained stable (Fig. 3A), and substrate type, water velocity, and depth determined the suitability scores (Fig. 4A). For this particular case, low-flow conditions were within the range of tolerance described in the habitat suitability rules and most of the domain appeared to be appropriate for mussels. High flows were the limiting conditions and substrate stability was the determining factor for medium to high flows (2039 m3/s and 3965 m3/s, respectively; Fig. 4B, C). The habitat suitability distribution at 3965 m3/s was used for subsequent delineation of suitable areas because it represents the expected yearly peak discharge and the 3.4-y return-period event, which is within the range suggested by Strayer (1999) as most sensitive for survival of freshwater mussels.
Formation of mussel beds
Flow patterns determined those areas where juveniles accumulated. Juveniles were deposited primarily in lower-velocity regions near the edge of the river and along the border of the main navigation channel (Fig. 5A). Higher velocity areas along the main and secondary channels generally prevented settlement. Those juveniles settling in suitable habitats survived and ultimately created mussel beds, but the others died because of inappropriate habitat conditions.
The simulation results indicate that juvenile mussels may travel from a few meters to several kilometers with the flow before settling (Fig. 5B). At a settling velocity of 0.03 cm/s, juveniles required >3000 s to drop 1 m in the water column. If they entered the fast-flowing regions in the main channel while dropping, they were transported at velocities on the order of 0.7 m/s (for an average flow of 2039 m3/s). At this velocity, they would travel ˜2 km in horizontal distance while dropping 1 m in the water column. Water depths in the main channel vary between 3 and 12 m. Therefore, many of the juveniles produced in the domain are likely to be transported downstream.
The simulated spatial distribution of mussel beds in this reach of the UMR (Fig. 5C) strongly coincided with reported locations of diverse and abundant mussel beds (USACE 1981, 1984). The only exception was a mussel bed located close to the upstream domain boundary and whose recruitment probably comes from upstream of the study area. The mussel locations are from the 1980s, but 3 of the beds were verified in 2004. Various mussel species were present (Amblema plicata being the most common), and the presence of juveniles indicated active recruitment. After 70 simulated years, estimated mussel densities varied between 5 and 98 mussels/m2, densities within the range of abundant mussel beds surveyed in a reach of the UMR 40 km upstream of the study area (Whitney et al. 1997). No records of mussel densities in the study area were found for comparison.
A numerical model was used to simulate the spatial distribution of developing mussel beds in a reach of the Mississippi River given a variety of hydrodynamic and substrate input parameters. The spatial distribution of the emerging mussel beds was strongly correlated with the known locations of mussel beds in this reach, suggesting that this model can accurately capture the relevant physical and biological processes driving population dynamics in a large river. To our knowledge, this represents the first such attempt with unionids.
The successful identification of mussel beds in the UMR was possible largely because of the application of the RSS. Several authors have predicted the distribution of unionids in small- to medium-sized rivers using various derivatives of shear stress (Layzer and Madison 1995, Hardison and Layzer 2001, Howard and Cuffey 2003). Unfortunately, the applicability of these results to other case studies is limited because typical shear stress measures are flow-conditional and the effect of shear force depends on substrate composition. Normalizing shear stress with the nondimensional parameter RSS and applying the rule RSS < 1 for substrate stability defined a measure that is not flow conditional and may be used across varying flow regimes and sediment types. Our data suggest that this concept is applicable in systems ranging from a 4th-order stream (Layzer and Madison 1995) to large rivers (our study). In theory, data from studies in rivers with different substrate and hydrodynamic characteristics could be normalized with RSS and compared or used subsequently to assess species-specific tolerances to substrate stability.
The computation of RSS at different flow rates showed how sediment stability and the location of flow refuges varied with flow rate, as previously noted by Strayer (1999). Large areas in the domain seem stable under low flows when, in fact, they experience active sediment motion that would render them unsuitable for mussels to inhabit at medium to high flows. During floods, the river bottom often shifts, and mussels may be crushed, buried, or sent downstream (Strayer 1999, Hastie et al. 2001). Periodic scouring events may determine the age structure of mussel populations in many rivers (Vannote and Minshall 1982). Nevertheless, we hypothesize that it is the annual peak flows that most often limit the spatial distribution of freshwater mussel communities. In our analysis, the computation of RSS at the mean annual peak flow provided an appropriate estimate of the location of flow refuges and the spatial distribution of suitable habitats for mussel survival. The RSS also may be used to predict areas where juvenile mussels cannot settle. Jowett (2003) commented that hydrodynamic conditions required for erosion of very fine particles can be considerably higher than those required for deposition. Juvenile mussels are small particles (on the order of 0.2 mm), so they are unlikely to be able to settle in areas where particles ≥0.25 mm are actively transported with the flow.
The RSS has application beyond the mussel dynamics simulations. It is analogous to the hypothetical habitat preferences for benthic invertebrates in gravel-bed rivers presented by Jowett (2003) that link substrate stability with mean velocity and substrate particle size. Nevertheless, RSS has 2 significant advantages over the Jowett (2003) model: 1) it is nondimensional and not flow-conditional, and 2) it provides a direct measure of substrate stability by addressing the cause (shear stress) and effect (sediment motion) relationship. As such, RSS may be applied in the evaluation of best management practices involving changes in flow patterns across a variety of river systems.
Simulation of juvenile dispersal
The simulation of dispersal of juveniles with flow was critical for determining which suitable habitats could be colonized by mussels and why areas with similar habitat suitability might have different mussel densities. The simulation of linear distances travelled by juveniles after they detach from their fish host and the spatial arrangement of mussel beds in this river reach suggest that mussel beds may be connected by recruitment events, at least under certain flow regimes. These data support the hypothesis of Vaughn and Taylor (2000) suggesting that initial colonization of habitats by unionids depends primarily on regional processes such as the dispersal of host fishes (i.e., dispersal of parasitic mussels), whereas growth and reproduction are most likely influenced by local environmental processes. Similarly, Vannote and Minshall (1982) hypothesized that block boulders prevented significant bed scour during major floods and these boulder-sheltered beds, although rare, contributed to population recruitment elsewhere in the river. The simulation results help us to understand why flow rates during glochidial release and juvenile settlement may limit juvenile recruitment (Hardison and Layzer 2001, Hastie et al. 2001) and, depending on the magnitude and duration, also may have long-term repercussions on population dynamics (Vannote and Minshall 1982, Di Maio and Corkum 1995). If this relationship between recruitment and flow rates holds true across a variety of systems, changes in hydrodynamic patterns may alter the evolution of mussel beds. Modifications of flow patterns that increase flow velocities may preclude recruitment of young individuals and hinder the long-term survival of otherwise healthy mussel beds.
The analysis presented here, while useful for predicting the evolution of mussel beds, is limited by resolution and quality of the input data. For instance, suitable patches smaller than the average cell size could not be captured by the mussel dynamics model. Nevertheless, given the scale of analysis, the results are reasonably accurate for a first attempt at modelling unionid populations in a large river by integrating ecological and hydrodynamic information. The mussel dynamics model can be a useful tool for evaluating the potential effects of different management practices on the overall evolution of freshwater mussel populations and can help to discriminate those areas with the highest potential to sustain healthy mussel communities and promote a more efficient allocation of resources to protect these areas.
This work was funded by IIHR–Hydroscience and Engineering, University of Iowa. The authors acknowledge Nathan Young and Andrew McCoy for conducting the fieldwork and developing the hydrodynamic model, Jeffrey Steuer for insightful comments on the applicability of the shear stress ratio concept, and WL Delft Hydraulics, the US Geological Survey, US Army Corps of Engineers, and Iowa Department of Natural Resources for their willingness to share data and provide guidance during the model development and validation.
Appendix. Preliminary verification of the shear stress ratio concept
Layzer and Madison (1995) investigated microhabitat use by unionids in a small stream in Kentucky, USA. The shear stress ratio (RSS) was applied to their data, and water depth and substrate preferences for unionids could be explained in terms of substrate stability (Fig. A). Silt substrates reached the onset of sediment motion at a RSS = 1 even at shallow water depths; hence, silt substrates were used only by a very small fraction of the mussels. Under low discharge conditions, mussels preferred depths between 0.1 and 0.4 m, the values for which RSS < 1 for all substrate classes, excluding silt. Under high discharge condition, water depths varied between 0.4 and 1.5 m. Mussels did not have access to low RSS areas, and their preference for water depths ranged from 0.75 to 1.0 m, corresponding to areas with RSS ≤ 1. RSS was consistently <1 for gravel substrates across all water depths. As a result, the relative number of mussels using gravel was disproportionately greater than for any other substrate type. Relative use of cobbles and boulders varied greatly among species because mussels also require small particles for burrowing (Layzer and Madison 1995).