Effects of vertical mixing on the Lake Michigan food web: an application of a linked end-to-end earth system model framework

Physical processes may affect ecosystem structure and function through the accumulation, transport, and dispersal of organic and inorganic materials, nutrients, and organisms; they structure physical habitat and can influence predator–prey interactions and trophic production. In the Laurentian Great Lakes, horizontal currents generally dominate, but little is known about the effects of vertical mixing on lake food webs. We developed a linked earth system model and used it to explore how vertical mixing affects the productivity of Lake Michigan (LM), the world’s fifth-largest lake, whose food web and fisheries have been adversely affected by invasive Dreissena mussels. We hypothesized that higher vertical mixing would result in higher food web biomass by making phosphorus more available to the lower food web, and that filtration by invasive mussels would counter the effects of mixing and decrease food web biomass. Using linked climate, hydrodynamics, and ecosystem models, we projected the response of LM’s food web to scenarios of different levels of vertical mixing, with and without invasive mussels. Biomass of most functional food web groups increased with increases in vertical mixing, with the greatest increases in phytoplankton and zooplankton. Increased biomass was due to the replenishment of nutrients into the euphotic zone, which enhanced growth and biomass of lower trophic levels through bottom-up effects. However, filtration by invasive mussels reduced the positive effects of mixing for most species. Future applications of the linked earth system framework will explore the effects of climate warming and nutrient reduction on fisheries production to inform fisheries managers.

found on the bottom into the water column. In temperate regions, the degree and effectiveness of mixing has a seasonal dependence; during the thermally stratified seasons (summer), the exchange between surface waters and deeper depths (hypolimnion) is restricted, but during the non-stratified seasons (late fall, winter, and early spring), the thermal conditions are isothermal, the water column is well mixed, and vertical exchange is high (Wetzel 2001). The importance of these seasonal mixing patterns has been demonstrated in one-and two-dimensional water quality model studies (e.g., Boegman et al. 2008); with increased computing power, 3-dimensional models with fine spatial resolution are now prevalent (e.g., Rowe et al. 2015Rowe et al. , 2017Chen et al. 2002), which can emphasize both the horizontal and vertical physical processes that help to structure ecosystems.
The importance of VTD for ecosystem dynamics in the Laurentian Great Lakes is further magnified by its effects on the food web that are modified by the invasive dreissenid mussels (quagga mussel Dreissena rostriforma bugensis, and zebra mussel D. polymorpha). Dreissenid mussels are benthic filter-feeders that have altered water quality, habitat, biogeochemistry, and energy flow in four of the five Great Lakes and have negatively affected food webs and fisheries (Kerfoot et al. 2010;Vanderploeg et al. 2010;Madenjian et al. 2015;Boucher 2019). The growth of dreissenid mussels depends on the availability of organic material obtained through passive sinking from the pelagic zone, or more often, delivery through VTD to the bottom. Moreover, dreissenid mussels affect the cycling and offshore transport of phosphorus (Hecky et al. 2004). Despite phosphorus loadings to Lake Michigan remaining relatively constant, water column concentrations of nutrients have declined after quagga mussels irrupted in LM in 2004, particularly in offshore waters. Dreissenid mussel filtration acted to sequester and divert nutrients from pelagic waters to the bottom and nearshore (Hecky et al. 2004), thereby depleting nutrient concentrations available to the offshore pelagic food web. Thus, how vertical and horizontal transport processes affect dreissenid mussel filtration, growth, and food web effects are of interest to lake stakeholders, resource managers, and limnologists (Rowe et al. 2015;Shen et al. 2018Shen et al. , 2020 and necessitate the incorporation of VTD in a 3-dimensional ecosystem model. Clearly, hydrodynamics plays a major role in many ecosystem processes, but they are often poorly parameterized in ecosystem models, especially those which include dynamics of an upper food web. For example, in the Atlantis Ecosystem Model (Fulton et al. 2004a(Fulton et al. , 2004b, a popular trophodynamics ecosystem model framework, VTD is largely ignored (Audzijonyte et al. 2017), yet it is an important and potentially dominant component of vertical transport materials and nutrients in large lakes and estuaries, with mixing timescales that are several orders of magnitude smaller than those associated with advection or settling (e.g., Koseff et al. 1993). While transient wind-induced upwellings occur in many large lakes, including the Laurentian Great Lakes (e.g., Plattner et al. 2006;Rao and Schwab 2007;Li et al. 2021), VTD likely plays a more dominant role in vertical transport over large temporal and spatial scales, especially during periods of weak and no thermal stratification which can occur for up to two-thirds of the year. Thus, accurate representations of hydrodynamic processes, including seasonal mixing cycles, are imperative for meaningful ecosystem modeling and forecasting.
Herein, we evaluate the response of the Lake Michigan food web to VTD. Specifically, we focused on four VTD treatments and their influence on food web effects with and without invasive dreissenid mussels. To accomplish this, we used our Great Lakes Earth System Model (GLESM), which is a linked model framework that uses outputs of downscaled temperature, wind, and precipitation from a climate model as external forcings for a hydrodynamics model (Great Lakes Finite Volume Community Ocean Model, GL-FVCOM), then takes the hydrodynamic output and nutrient estimates from major tributaries (Luo et al. 2017) to drive an end-toend ecosystem model-the Lake Michigan Atlantis Model (LMAM).

Study ecosystem
The North American Laurentian Great Lakes (Fig. 1) are the largest freshwater system on the planet, with a total surface area of 244,160 km 2 , and contain about 21% of the world's surface freshwater and about 84% of North America's surface freshwater (https:// ijc. org/ sites/ defau lt/ files/ SAB_ Great-Lakes-Scien ce-Strat egy-summa ry-report_ 2022. pdf accessed on 1/3/2023). The lakes also provide valuable ecosystem services to over 38 million people residing in the basin. The lakes provide drinking water and support agriculture, transportation, commercial shipping, hydroelectric power, and a wide variety of recreational opportunities including a world-class fishery that generates approximately $7 billion USD per year (Krantzberg and De Boer 2008). Of the five Great Lakes, Lake Michigan is the second largest Great Lake by volume (4918 km 3 ), the third largest by area (57,763 km 2 ), and is also the fifth largest lake (by area, 6th by volume) in the world (https:// www. glc. org/ lakes/ lake-michi gan, accessed 12/21/2022). As with many large ecosystems, Lake Michigan has, and continues to experience, many humaninduced stressors including climate change, invasive species, and high nutrient concentrations in coastal waters that contribute to harmful algal blooms (HABs) and hypoxia Tomlinson et al. 2010;Allan et al. 1 3 2013; Klump et al. 2018;Miller et al. 2022), but low nutrient concentrations in the offshore area.

Great Lakes Earth System Model (GLESM): linked model framework overview
GLESM integrates and links climate, weather, watersheds, hydrodynamics, and chemical and food web processes into a modeling framework that provides capabilities for hindcasts, nowcasts, and forecasts of regional effects (natural and anthropogenic) on the Laurentian Great Lakes. Here, we highlight GLESM's capabilities through weather-forced circulation and mixing on the food web of the Laurentian Great Lakes.

Hydrodynamics
Lake Michigan's physics and thermodynamics were simulated using the Great Lakes Finite Volume Community Ocean Model (GL-FVCOM), a spherical-coordinate, unstructured-grid hydrodynamic model designed to solve the primitive (i.e., hydrostatic) equations. GL-FVCOM was originally implemented by Bai et al. (2013), who modified the original FVCOM initially developed by Chen et al. (2003Chen et al. ( , 2006Chen et al. ( , 2013. GL-FVCOM incorporates several modifications to the original FVCOM scripts (version 3.1.6), including the use of the internal ice module (Gao et al. 2011) based on a coupled elastic-viscous-plastic rheology ice model (CICE: Hunke and Dukowicz 1997) and a wind-wave mixing boundary condition scheme developed by Hu and Wang (2010). Both modifications dramatically improved modeled surface temperature and subsurface thermal structure in the Laurentian Great Lakes, which are particularly important for ecological modeling. The horizontal computational grid is composed of unstructured triangular elements, and vertical layers are delineated using a terrain-following σ-coordinate scheme with increased resolution at the surface and bottom boundaries. Additional details on model parameterization can be found in Cannon et al. (2023) and Wang et al. (2023). For this study, the GL-FVCOM model domain included both Lake Michigan and Lake Huron, which are connected through the Straits of Mackinac and share a common water elevation (Fig. 1). Although the lakes are often discussed and managed separately, their hydrologic connection requires combined modeling for accurate hydrodynamic simulations. For simplicity, the lakes were modeled as an enclosed basin, and all mass fluxes (i.e., inflows, outflows, evaporation, precipitation) were ignored. The horizontal resolution of the model grid ranged from 3 to 5 km, with 21 vertical sigma levels used for subsurface computations. Model results for each lake were separated during post-processing, and only model grid cells in Lake Michigan were retained for further analysis. Model outputs, including vertical and lateral advective exchange rates and water column temperatures, were rescaled and averaged over larger sub-domains and depth layers (see Supplementary Information) for input into the Atlantis Ecosystem Model (discussed below: 34 boxes; maximum 6 layers per box).

Ecosystem model
Lake Michigan Atlantis Model (LMAM) was developed and integrates physical, chemical, ecological, and fisheries dynamics in a three-dimensional, spatially explicit framework. The framework was initially developed by Fulton et al. (Fulton 2001;Fulton et al. 2003Fulton et al. , 2004aFulton et al. , 2004bFulton et al. , 2004c) that includes a 3-dimensional rendering of an ecosystem and five modules: hydrographic, ecology, and three management and human use modules (harvest sub-model, assessment submodel, and economics sub-model). Here, we only focus on the two core modules: hydrographic and ecology modules.
The hydrographic module in LMAM accepts inputs (i.e., water temperature, vertical and horizontal fluxes) from the . GL-FVCOM simulated water temperatures are rescaled to the average for each polygon and depth stratum, while fluxes (vertical and horizontal) are the average flux (direction and magnitude) across each polygon/stratum segment rescaled output of GL-FVCOM, which are used as forcing functions in the ecological module. In the ecological module, nutrient dynamics, organismal growth, mortality, consumption, recruitment, physical transport/behavioral movement, and habitat dependency are modeled in each polygon as a system of differential equations on a 12-h time step. Primary producers and invertebrates are represented as biomass pools and vertebrates as age-structured populations (Table 1).
Here, our emphasis is on vertical mixing and its effects on the food web. Thus, modifications to vertical mixing were done in the hydrographic module of Atlantis and not in GL-FVCOM. Atlantis includes a simple diffusion function designed to simulate mixing throughout the water column using a constant turbulent vertical diffusion coefficient (Kz), which is often set to molecular or nearmolecular diffusion levels (1 × 10 −7 − 10 −6 m 2 /s). These weak mixing coefficients are unrealistic in Lake Michigan, where vertically averaged hypolimnetic mixing rates have been observed to vary between 1 × 10 −5 and 10 −2 m 2 /s over the seasonal stratification cycle (Cannon et al. 2021). To explore this assumption of weak vertical mixing, we designed four parameterized mixing treatments: (1) three levels of annual average constant VTD (Kz = 1 × 10 −5 , 10 −4 , 10 −3 m 2 /s; hereafter referred to as C5, C4, and C3, respectively); and (2) a sinusoidally varying seasonal VTD (hereafter referred to as S, with an annual average of 1 × 10 −3.85 m 2 /s) ( Table 2, Fig. 2). The S scenario was designed to mimic the modeled (and observed) seasonal cycle of mixing in Lake Michigan (Fig. 2), with depthaveraged turbulent diffusion rates that varied between 1 × 10 −5 m 2 /s in the summer and 1 × 10 −2 m 2 /s in the winter. S coefficients (K z-seasonal ) were explicitly defined based on the day of year (doy, range: 1-366), such that Table 1 List of model groups in the Lake Michigan Atlantis Model. Model groups are categorized as "functional groups," including piscivorous fish (PISC), omnivorous fish (OMNI), planktivorous fish (PLAN), benthivorous fish (BENV), benthic invertebrates (BENT), zooplankton (ZOOP), macro-algae (MACR), phytoplankton (PHYTO), microbes (MICR), detritus (DETR). Initial biomass (MT/ km 2 ) is the basin wide average biomass at the start of the model simulation based on a Lake Michigan Ecopath model (Rutherford et al. 2021). "Modeled as" refers to whether the model group was modeled as an age-structured population or as a biomass pool For each mixing scenario, VTD coefficients in Atlantis were set equal to Kz = 1 × 10 −5 , 10 −4 , 10 −3 m 2 /s or K z-seasonal at all depth layers in offshore model boxes (i.e., boxes not touching the coast). For all scenarios, nearshore (i.e., boxes along the coast) VTD coefficients were increased by an order of magnitude to simulate the turbulence-enhancing effects of coastal boundary layers.
We recognize that because the seasonal VTD coefficient described here is not directly linked to hydrodynamic forcing inputs (i.e., lake temperatures, current velocities), the day-to-day lake physics may be dynamically inconsistent. For example, there may be days where parameterized mixing rates are high, despite the presence of turbulence-limiting thermal stratification in a given model box. However, given the relatively consistent annual cycle of lake-wide vertical mixing rates (e.g., Fig. 2), we expect the effects of dynamic inconsistencies to be relatively minor, especially when integrated over the large spatial and temporal scales discussed in this manuscript. (1)

S_M*
Sinusoidal seasonal vertical mixing, mussels S_noM Sinusoidal seasonal vertical mixing, no mussels C5_M Constant vertical mixing (K z = 1 × 10 −5 m 2 /s), mussels C5_noM Constant vertical mixing (K z = 1 × 10 −5 m 2 /s), no mussels C4_M Constant vertical mixing (K z = 1 × 10 −4 m 2 /s) mussels C4_noM Constant vertical mixing (K z = 1 × 10 −4 m 2 /s), no mussels C3_M Constant vertical mixing (K z = 1 × 10 −3 m 2 /s), mussels C3_noM Constant vertical mixing (K z = 1 × 10 −3 m 2 /s), no mussels The S parameterization was for seasonally varying VTD, while C3-C5 were constant VTD parameterizations. The C5 VTD (red) was the baseline for comparison of VTD effects lake-averaged biomass of 21 model groups or species. We focused on the calibration of phytoplankton, the base of the food web, and on upper trophic levels, especially alewife, chinook salmon, coho salmon, and lake trout. These groups are key drivers of Lake Michigan's economically valuable recreational fisheries and a current focus of fisheries managers and contain the most complete time series data.

Scenario simulations and analysis
We used LMAM to run eight scenario simulations to understand how VTD may affect the food web biomass, with and without the presence of invasive mussels in the food web (Table 2). Specifically, the scenarios included eight full combinations of two mussel treatments (with mussels and without mussels) and four VTD treatments, and each scenario was run with a 12-h time step for 100 simulation years. We compared the model output of the vertical distribution of SRP and nanoplankton in a deep offshore habitat, box 23 ( Fig. 1) during winter and summer. We also calculated the total phosphorus in the water column of box 23 and in the euphotic zone (to 40 m) for each scenario during winter and summer. Furthermore, we aggregated the biomass of model species or groups within similar trophic levels together into 9 functional groups, including piscivorous fish, omnivorous fish, planktivorous fish, benthivorous fish, benthic invertebrates, zooplankton, macroalgae, phytoplankton, and microbes (Table 1). For reporting purposes, simulated biomasses of functional groups were averaged across horizontal polygons and depths for the whole lake (MT/km 2 ) and summarized as average biomass for the last 20 years (years 80-100) of the simulation. The biomass response of functional groups to the VTD mixing scenarios was compared to functional group biomass response to the corresponding lowest VTD mixing scenarios (C5_M or C5_noM as baselines) and was reported as the relative percent change as follows, using seasonally varying mixing with mussels S_M as an example): We considered biomass changes of functional groups relative to the corresponding baseline C5 scenarios (C5_M and C5_noM) as strong effects if their absolute values were ≧ 20%, as intermediate effects if absolute change values were from 10 to 19%, and as negligible effects if absolute change values were < 10%.

Model calibrations to data (GL-FVCOM, Atlantis)
GL-FVCOM simulated forcing data were validated against observations in Lake Michigan. Lake-averaged ice cover and lake surface temperatures were compared to remote sensing observations available through the Great Lakes Surface Environmental Analysis (https:// coast watch. glerl. noaa. gov). Ice cover estimates were generally within 10% of observations, and mean lake surface temperature biases were less than 0.75 ℃, with a root mean square error of 1 ℃ in Lake Michigan. Modeled temperature structure and stratification were compared to in situ mooring observations in each lake (NOAA GLERL 2019a, 2019b). Modeled temperature structure highlighted overly diffuse thermoclines during the stratified period, with maximum root mean square errors of ~ 3 ℃ at 30 m depth. These biases are not expected to have a detrimental impact on ecological simulations, especially considering the relatively low vertical resolution of the ecological model (6 layers) compared to the hydrodynamic model (21 layers). Simulated seasonal surface current patterns (Figure SI-1) also compared favorably with previous observational (Beletsky et al. 1999) and modeling (e.g., Beletsky and Schwab 2008) studies, with increased coastal current speeds and counter-clockwise gyres in the north and south basins. Additional details on model validation can be found in Cannon et al. (2023). Atlantis model calibration was performed for 21 food web groups ( Figures SI-6

VTD effects on vertical distribution of nutrients and nanoplankton
Nutrients and plankton showed different vertical distributions under VTD scenarios (Fig. 3). Phosphorus is often the limiting nutrient in freshwater, so we present a vertical distribution of soluble reactive phosphorus (SRP) in a representative offshore model box (box 23, located in southeast Lake Michigan). For scenarios of C5_M, C4_M, C3_M, and S_M, SRP concentrations during summer were low in the upper water column and high in the lower water column, with a very low SRP concentration at the surface. SRP concentration in the upper water column was the highest under the C3_M scenario. During winter, SRP concentration increased in the water column compared to summer and was distributed almost homogeneously in the water column under S_M and C3_M. However, for C5_M and C4_M during the winter, SRP showed a similar pattern as the summer, lower in the upper water column and higher in the lower water column. Nanoplankton biomass exhibited a different seasonal vertical distribution than did SRP. During summer, nanoplankton biomass was high in the upper water column and low in the lower water column under the scenarios of S_M and C3_M, but was consistently low throughout the water column under the scenarios of C5_M and C4_M. During winter, nanoplankton biomass was homogenous and low throughout the water column under all four mixing scenarios.
The total phosphorus (TP) content in the water column (SRP plus phosphorus from organic sources including organismal biomass) increased with increased levels of constant mixing (from low C5_M to high at C3_M), while S_M had a similar water column TP concentration as under C4_M for summer and winter. TP concentration in the water column was similar or slightly higher during summer than during winter for each scenario (Table 3). TP concentration in the euphotic zone (< 40 m) also increased with increased levels of constant mixing and was relatively stable between summer and winter across scenarios.

VTD effects on the Lake Michigan food web with invasive mussels
Under the treatments with invasive mussels, compared to the low constant mixing scenario (C5_M), the total biomass of the whole food web increased by 24, 42%, and 541% with increasing levels of mixing (under 4C_M, S_M, and 3C_M, respectively) ( Table 4). When functional food web groups were considered, the C4_M mixing scenario had large positive effects (≥ 20% above the baseline C5_M biomass) on the  biomass of omnivorous fish, zooplankton, macroalgae, and phytoplankton, increasing by 85, 64, 26, and 61%, respectively, and had intermediate positive effects (10-19% above C5_M biomass) on piscivorous fish biomass by 11%, but had negligible changes (between − 10% and + 10% change from C5_M biomass) on planktivorous fish, benthivorous fish, benthos, and microbes. The S_M mixing scenario had similar positive effects on the functional groups as C4_M, except that S_M increased microbe biomass by 20%. Both C4_M and S_M showed negligible negative effects on benthivorous fish and benthos. The highest mixing scenario (C3_M) had large positive effects on all functional groups, except on benthivorous fish which increased by 11% above baseline.

VTD effects on the food web without dreissenid mussels
Under the no-mussel treatment, compared to baseline C5_ noM, the total biomass of the food web under the C4_noM, S_noM, and C3_noM scenarios increased by 85, 97, and 734%, respectively (Table 4). For the aggregated trophic groups, C4_noM, S_noM, and C3_noM had large positive effects on almost all functional groups, except that C4_noM had negligible effects on planktivorous fish and benthos, and S_M had intermediate positive effects on planktivorous fish and benthos. Compared to other no-mussel scenarios, the C3_noM scenario had the largest positive effects on each functional group.

Food web effects of VTD between scenarios of invasive mussels and no mussels
Relative to the baseline C5 scenarios, higher mixing scenarios increased the total biomass of the food web more under the no-mussel treatment compared to a treatment with mussels in the food web (85% vs 24%, 97% vs 42%, and 734% vs 542% for C4, S, and C3, respectively). For all functional groups, percent changes in biomass were higher for the scenario of no invasive mussels in the food web compared to one with invasive mussels in the food web (Table 4). Moreover, C4 and S mixing scenarios with mussels had slightly negative effects (although negligible) on benthivorous fish and benthos, but had positive effects (up to 62% increase in biomass) under the scenario of no invasive mussels.

Discussion
Our scenario simulations indicated that VTD mixing had an overall positive effect on most food web groups, with some exceptions. VTD acted to replenish nutrients in the euphotic zone which stimulated production up the food web from phytoplankton to microbes, then to zooplankton and their invertebrate and vertebrate predators, and ultimately to piscivores. Differences in magnitude and direction of organism biomass response to VTD were also driven by invasive mussel filtration of phytoplankton and alteration of nutrient cycling, and predator prey interactions which acted to dampen the biomass response of some species or groups to VTD. Effects of VTD mixing generally were greater for lower trophic levels than for higher trophic levels and greater for pelagic organisms than for benthic organisms. Most benthic invertebrates and their predators, benthivorous fish had a negligible response to low VTD scenarios. The pelagic food web response to VTD was higher for phytoplankton and zooplankton than for planktivorous fish and their predators, piscivorous fish (Leonhardt et al. 2020). The difference between upper and lower trophic levels in their response to VTD and mussels is consistent with observations by McMeans et al. (2016) who suggested that the generalized feeding strategies and movements of large, higher trophic taxa can lead to a food web structure that is adaptive in its response to environmental changes such as variability in production or effects of invasive species.
Several studies have documented negative effects of invasive mussels on food web biomass and trophic transfer efficiency (Bunnell et al. 2018;Trumpinkas et al. 2022). Our model results suggest that ecosystem biomass may be further suppressed under a climate warming scenario, which would weaken VTD by increasing the duration of thermal stratification and shift the lake from a dimictic condition to a more monomictic condition. Rowe et al. (2017) modeled the relative influence of climate warming, nutrients, and invasive mussels on the Lake Michigan lower food web and found that the influence of Dreissena mussels on plankton was reduced during thermal stratification. Their model suggested warmer weather would deepen mixing of the water column, thereby reducing nutrient flux and increasing mussel filtration in the nearshore but lower mussel filtration offshore.
To build a framework of complex models to address critical issues of environmental change, we made assumptions and accepted variance in model tuning that we plan to address in future applications. The Atlantis model calibration was satisfactory for many but not all food web species or groups. For example, we may have underestimated the effects of invasive mussels on the food web response to VTD. The LMAM coupled framework was able to project the relative influence of VTD on invasive mussel cycling of phosphorus and reduction in biomass of phytoplankton caused by dressenids grazing in Lake Michigan food web, which has been well documented for the Great Lakes food webs (Hecky et al. 2004;Fahnenstiel et al. 2010;Vanderploeg et al. 2015;Pothoven and Vanderploeg 2020). However, not included in our simulations were consequences of the indirect effects of mussels on light, the vulnerability of prey to predators, consequent behavioral responses and interactions (Bourdeau et al. 2015), or the increased importance of mussel veligers to plankton, larval fish diets, and fish recruitment in the food web (Eppehimer et al. 2019).
Our study demonstrated the effects of vertical mixing on the Lake Michigan food web could be similar between using a model/observation-based sinusoidally varying seasonal VTD or an annual averaged constant VTD. Using an annual averaged constant VTD may be reasonable for investigating annual phenomena, but a model/observationbased sinusoidally varying seasonal VTD has its advantages in projecting seasonal phenomena or those of finer temporal scales, and these seasonal measures can provide the means to help refine the GL-FVCOM model for determining longer-term phenomena associated with climate change. Moreover, seasonal varying VTD is how a temperate lake operates. In the Great Lakes, VTD varies seasonally (e.g., Cannon et al. 2021) with limited diffusion rates (Kz = 1 × 10 −5 m 2 /s) during the thermally stratified summer months and with energetic and full water column mixing rates (Kz = 1 × 10 −4 − 10 −2 m 2 /s) during the fallwinter-early spring isothermal months.
We also exaggerated the difference among mixing treatments in the nearshore zone. When lakes are covered by ice, there should be less mixing in winter, but still full water column mixing in the spring and fall. Ice cover can drastically reduce vertical mixing during winter, blocking wind forcing and generating turbulence-suppressing inverse stratification in the surface mixed layer. As well, there have been few observations of Great Lakes ecosystems taken under the ice, and there are virtually no data on how the Great Lakes ecosystems respond to vertical mixing under winter conditions Ozersky et al. 2021). In future, we plan to examine the effects of reduced ice cover on vertical mixing to understand how that may affect fish species recruitment and food web productivity.
Although it was beyond the scope of the present study, which focused on testing the effects of improved turbulence parameterizations, in the future, we plan to modify the LMAM to accept time series of vertically variable turbulent diffusion coefficients as forcing inputs from more observations or from hydrodynamics outputs. While assumptions of depth-uniform vertical turbulent diffusion coefficients are reasonable for the energetic fall and winter, the vertical structure of VTD is heavily influenced by water density stratification during the summer, when mixing rates outside the surface and bottom boundary layers fall to near-molecular levels (i.e., K z ≈ 10 −7 m 2 /s; Cannon et al. 2021). This weak hypolimnetic mixing acts as a barrier to benthic particle delivery, reducing the effective filtration rate of mussels in deeper waters (Shen et al. 2018;Rowe et al. 2017).

Conclusions
1. We used linked climate, hydrodynamics, and ecosystem models to project the response of LM's food web to seasonally variable and seasonally constant vertical mixing scenarios, with and without invasive mussels. 2. Relative to constant low vertical mixing, seasonally varying vertical mixing increased biomass of the total food web biomass by 42% and of most food web groups by replenishing nutrients in the euphotic zone, which enhanced phytoplankton growth and biomass of lower trophic levels through bottom-up effects. High levels of constant vertical mixing acted to further increase the biomass of most food groups. 3. Without invasive mussels in the food web, all scenarios of vertical mixing acted to increase the biomass of the total food web biomass. Filtration by invasive mussels reduced the positive effects of seasonally varying mixing on most species. 4. A constant vertical mixing scenarios of Kz = 1 × 10 −4 m 2 /s provided similar food web effects as the seasonally variable mixing scenario. However, (a) seasonally varying vertical mixing best represents what actually occurs in nature (Cannon et al. 2021), and it can be captured in hydrodynamic models, and (b) the ability of a single constant vertical mixing coefficient may vary annually, and the ability of the model to capture future conditions under climate change cannot be assumed. Therefore, we recommend using seasonal varying vertical mixing to best capture current and future food web response to seasonal and annual changes in lake hydrodynamics, water temperature, precipitation, and ice cover resulting from climate change.