For lakes experiencing extended ice-cover seasons, ice phenology has a substantive impact on thermal structure and dissolved oxygen (DO) dynamics during winters. This study applied a three-dimensional lake model (AEM3D) in Lake Winnipeg over 2016-2018, spanning two full winter seasons. Sensitivity analysis showed that the modeled ice cover thickness, formation, and duration were most sensitive to snow depth and snow/ice albedo. The model well simulated the ice freeze-up timing with less than 5 days discrepancy, but it underestimated the ice cover thickness and overestimated the ice cover duration in the unusual warm winter (2016-2017). Inverse stratification was developed under the ice, but the model could not fully reproduce it due to a lack of a sediment heat flux component. DO decreased with the formation of ice cover, leading to bottom hypoxia over the lake. The model indicates that ice phenology (i.e., ice cover duration, and blue/white ice thickness) affects the extent of winter hypoxia. We observed the DO decreased to < 2 mg L-1 (i.e., anoxia) in the North Basin, along with inverse stratification forming near lakebed, and an oxygen depletion rate reached 0.14 mg L-1 d-1 in the winter of 2016-17. The model captured but underestimated DO decline near the lakebed and simulated around 8% and 70% of lake area reached anoxia in winters of 2016-17 and 2017-18, respectively. This work provides insight into ice formation, under-ice thermal structure, and winter oxygen concentrations in a large prairie lake, and the role of changing ice phenology on northern lake ecology.

Parameterizations for bottom shear stress are required to predict sediment resuspension from field observations and within numerical models that do not resolve flow within the viscous sublayer. This study assessed three observation-based bottom shear stress (τb) parameterizations, including (1) the sum of surface wave stress and mean current (quadratic) stress (τb= τw +τc); (2) the log-law (τb= τL); and (3) the turbulent kinetic energy (τb= τTKE); using two years of observations from a large shallow lake. For this system, the parameterization τb= τw +τc was sufficient to qualitatively predict resuspension, since bottom currents and surface wave orbitals were the two major processes found to resuspend bottom sediments. However, the τL and τTKE parameterizations also captured the development of a nepheloid layer within the hypolimnion associated with high-frequency internal waves. Reynolds-averaged Navier-Stokes (RANS) equation models parameterize τb as the summation of modeled current-induced bottom stress (τc,m) and modelled surface wave-induced bottom stress (τw,m). The performance of different parameterizations for τc,m and τw,m in RANS models was assessed against the observations. The optimal parameterizations yielded root-mean-square errors of 0.031 and 0.025 Pa, respectively, when τc,m, and τw,m were set using a constant canonical drag coefficient. A RANS-based τL parameterization was developed; however, the grid-averaged modelled dissipation did not always match local observations, leading to O(10) errors in prediction of bottom stress. Turbulence-based parameterizations should be further developed for application to flows with mean shear-free boundary turbulence.