Environmental and climatic factors affecting winter hypoxia in a freshwater lake: evidence for a hypoxia refuge and for re-oxygenation prior to spring ice loss

Low dissolved oxygen, or hypoxia, is a common phenomenon in ice-covered lakes in winter. We measured dissolved oxygen (DO) before, during, and after ice-over to characterize the timing, severity, and spatial variability of winter hypoxia in Upper Red Rock Lake, Montana, home to one of the last remaining lacustrine populations of endemic Montana Arctic Grayling (Thymallus arcticus). Unlike most previous investigations of winterkill-prone lakes, we observed considerable horizontal spatial variability in DO, a non-linear winter oxygen depletion rate, and lake-wide re-oxygenation 2–4 weeks prior to spring ice loss. Parts of the upper 1 m of the lake and near stream mouths remained well-oxygenated even during late winter. DO levels were strongly associated with maximum daily air temperature. Our analysis of a 28-year weather record revealed large interannual variability in risk of winter hypoxia, with a slight declining trend in winter severity (number of days with maximum air temperatures ≤ 0°C) in Upper Red Rock Lake. The approach we used in our study provides a useful framework for quantifying and mapping the seasonal dynamics of the extent and severity of winter hypoxia, and for identifying critical winter habitats.


Introduction
Low dissolved oxygen (DO), or hypoxia, is a common and widespread phenomenon in ice-covered lakes in winter (Greenbank, 1945;Agbeti & Smol, 1995;Hasler et al., 2009). The onset of ice cover begins a period of physical separation of water and the atmosphere, which slows or halts processes that replenish DO. Ice and snow cover prevents oxygenation of lake water by diffusion and aeration, slows convective mixing, and reduces the amount of light available for photosynthetic production of oxygen by algae and submerged vegetation (Welch & Kalff, 1974;Malm et al., 1998;Terzhevik et al., 2009). If ice cover persists for a sufficient duration, the consumption of oxygen by microbial decomposition of dead and decaying organic matter and by absorption from bottom sediments may deplete DO to the point of hypoxia. In contrast, photosynthesis and DO production may continue to occur if light can still penetrate ice and snow cover. In most instances, the winter oxygen-limiting period ends with the breakup of lake ice in spring (Barica & Mathias, 1979).
Hypoxia can be a major cause of winter mortality in fishes (Greenbank, 1945;Cooper & Washburn, 1949;Hughes, 1973), which can limit species abundance and persistence (Danylchuk & Tonn, 2003) and fish assemblage composition (Tonn & Magnuson, 1982). The degree of winter hypoxia in ice-covered lakes is largely a function of pre-ice total oxygen storage, lake morphology, water exchange rate, duration of ice cover, extent of under-ice photosynthesis, and presence of inlet streams (Greenbank, 1945;Mathias & Barica, 1980;Guenther & Hubert, 1991;Leppi et al., 2016). Winter hypoxia is common in shallow, eutrophic lakes with high rates of winter oxygen depletion (Barica & Mathias, 1979;Mathias & Barica, 1980;Babin & Prepas, 1985;Terzhevik et al., 2009). Annual variability in winter air temperature and snowfall, however, may also have a significant effect on hypoxia development by influencing duration of the ice-cover period and depth of surface snow accumulation, which limits light penetration thereby reducing photosynthetic production of oxygen by algae and submerged aquatic vegetation (Barica & Mathias, 1979;Baehr & DeGrandpre, 2002).
Knowledge of the factors that contribute to the development of winter hypoxia allow for the identification of lakes at risk for winterkill and for the application of predictive models of winter DO concentration (Barica & Mathias, 1979;Fang & Stefan, 2000, White et al., 2008Leppi et al., 2016). The ability to predict whether a lake is at risk of severe winter hypoxia is an important management tool allowing preemptive action to avert large-scale winterkill of important fishery resources (Barica & Mathias, 1979;Barica et al., 1983;Prowse & Stephenson, 1986). Many predictive models of whole-lake DO conditions rely on lake morphometric parameters to predict DO concentrations and the likelihood of winterkill (Barica & Mathias, 1979;White et al., 2008;Leppi et al., 2016). Some winterkill risk models have incorporated weather conditions and landscape-level variables, in addition to the aforementioned lake attributes, to predict winter hypoxia severity (Fang & Stefan, 2000). To date, winter hypoxia models have generally relied on DO sampling from a single site to characterize whole-lake hypoxia severity in small lakes with limited surface inflow. However, the extent of winter hypoxia may be more heterogeneous, especially in larger lakes with surface inflow, requiring more extensive sampling to accurately characterize the spatial and temporal extent of hypoxia (Guenther & Hubert, 1991;Baehr & DeGrandpre, 2002;Hasler et al., 2009).
Our study focused on Upper Red Rock Lake, Montana, a shallow, eutrophic, high-elevation lake with high potential for development of hypoxia during winter ice cover (Gangloff, 1996). The lake is also home to Montana Arctic grayling (Thymallus arcticus), a now rare, glacial relict population confined to the Upper Missouri River watershed, and which occupy only a small fraction of their historic range (Liknes & Gould, 1987;Kaya, 1992). The population in Upper Red Rock Lake is the only significant remaining native, lacustrine population (Nelson, 1954;Kaya, 1992;Rens & Magee, 2007), and is genetically distinct from other Arctic grayling in Montana and from populations further north in Canada, and is in extremely low numbers (\ 300 spawners) (Peterson & Ardren, 2009;Warren et al., 2020). In a previous investigation, Gangloff (1996) observed widespread hypoxia in the lake, including evidence of fish mortality under the ice. Low overwinter survival due to winter hypoxia has been hypothesized as a limiting factor for this important population (Vincent, 1962;Gangloff, 1996;Mogen, 1996;Warren & Jaeger, 2017;Warren et al., 2020). However, Gangloff (1996) reported that hypoxia extent varies considerably due to annual fluctuations in environmental and weather conditions, and that inlet stream mouths may provide refugia from hypoxic conditions. However, the spatial extent and dynamics of winter hypoxia is not well known for this system. The objectives of this study were to assess the winter hypoxia threat in Upper Red Rock Lake by (1) determining the timing, severity, and spatial variability of winter hypoxia in Upper Red Rock Lake over two study years with varying weather conditions; (2) further refining the location and duration of potential refugia from hypoxia; (3) developing a predictive model of environmental factors affecting under-ice DO; and (4) examining historical trends in climate-related hypoxia risk.

Study site
Upper Red Rock Lake ( Fig. 1; 44°36 0 41 00 N, 111°43 0 36 00 ) lies at an elevation of 2017 m in the Centennial Valley in the Upper Missouri River Basin of southwest Montana. The 893 ha lake is surrounded by extensive wetlands and lies within Red Rock Lakes National Wildlife Refuge. It is a cold, polymictic lake that rarely stratifies during the spring-to-fall seasons due to shallow depth and frequent wind mixing (Cutting et al., 2016). The lake is fed by five inlet tributaries ( Fig. 1): two are reliant on snowmelt runoff from adjacent mountains (Red Rock and Tom creeks); two are spring-fed streams (Grayling, and East Shambow creeks); and one is a combination of upstream lake outlet and spring-fed discharge (Elk Springs Creek) (Gangloff, 1996). Flows exit the lake through a single outlet in the northwest corner of the lake that connects Upper Red Rock Lake to Lower Red Rock Lake and Red Rock River. Upper Red Rock Lake exhibits many characteristics typical of a winterkill-prone lake. The lake experiences long winters, typically remaining icecovered from November to April. Depth is shallow and rather uniform, with a maximum depth of 2 m during the ice-free period. Abundant rooted macrophytes (Potamogeton spp., Elodea spp., Myriophyllum spp.) exist on a uniform bottom comprised of mud, peat, and detritus (Gangloff, 1996;Mogen, 1996). The lake supports a diverse fish assemblage, including Arctic Grayling, burbot (Lota lota), white sucker (Catostomus commersonii), longnose sucker (Catostomus catostomus), longnose dace (Rhinichthys cataractae), mottled sculpin (Cottus bairdi), brook trout (Salvelinus fontinalis), and cutthroat-rainbow trout hybrids (Oncorhynchus clarkii bouvieri 9 Oncorhynchus mykiss; Gangloff, 1996;Cutting et al., 2016).

Lake sampling
To assess the extent, duration, and annual variability in lake DO during winter, measurements were made every 2-4 weeks during the winters of 2013-2014 and 2014-2015 from before ice-up (October) until travel on lake ice became unsafe during spring ice breakup. This resulted in 14 sampling occasions spread across the two winters: 8 in 2013-2014 and 6 in 2014-2015. Fifteen sample sites were generated randomly for each sampling occasion using the Create Random Points tool in ArcGIS (ver. 10.2; ESRI, Redlands, CA) and included both lake (10 sites) and stream mouth (5 sites) habitat types. Sample size was determined by the maximum number of sites that could be reliably sampled in a day, given frequent severe winter conditions that posed a safety risk to the sampling crew. Sample locations were stratified to detect anticipated differences between DO concentration in the lake and in areas \ 300 m of stream mouths (Gangloff, 1996). Holes were drilled into the ice with a 10-cm-diameter hand auger at each sample location, and DO (mg/l) and water temperature measured at 1-meter intervals from the ice-water interface to the bottom using an electronic DO probe (YSI Professional Plus, Yellow Spring, OH). Other possible covariates influencing under-ice DO that were also measured included ice thickness, snow depth, water depth, and air temperature. Air temperature data were obtained from a nearby weather station (Red Rock Montana NOAA weather station RRDM8; 44°40 0 59 00 N, 111°49 0 59 00 W, elevation 2,039 m) located about 12 km northwest and 22 m higher elevation than the lake.

Data analysis
At each sample site, DO and water temperature at 0, 1, and 2 m depth were averaged to generate mean values for each sample site. These values were then used to calculate lake-wide average DO concentration and water temperature by sampling occasion. Ice thickness and snow depth measurements from each sample site were used to calculate lake-wide averages in a similar manner. Variation in DO, water temperature, ice thickness, and snow depth across sampling occasions was assessed within years using a one-way ANOVA, followed by a Tukey-Kramer test to identify significant differences across sample dates. All statistical analyses were done with R programming language (R Core Development Team, 2018). The strength of evidence was classified as significant at the a \ 0.05 level.
To assess spatial variability in DO concentration in the lake, mean DO values for each sample site were calculated from measurements at 0 m and 1 m depth and input into ArcGIS. Only data from the upper 1 m of the water column were used as it was the only welloxygenated zone in the lake in both winters. The lower meter of the water column was consistently hypoxic (mean DO \ 2.0 mg/l in 2013-2014, \ 3.0 mg/l in 2014-2015) soon after complete ice cover. The upper meter, by contrast, was a comparatively dynamic layer, exhibiting a range of DO concentrations from * 0 to 10 mg/l. DO concentrations \ 4 mg/l were considered hypoxic in this situation because these DO concentrations in this range elicited stress responses (behavioral avoidance and increased opercular beats) and decreased survival in juvenile and adult Arctic Grayling in laboratory experiments whereas such impairments were absent at higher DO levels (Davis et al., 2019).
Maps of spatial distribution of DO concentration in this upper 1 m surface layer were determined for each sampling occasion using the inverse distance weighted interpolation tool in ArcMap Spatial Analyst based on interpolated raster datasets of 20 m 9 20 m grid cells of DO concentration. Accuracy of predicted DO from spatial mapping was assessed by comparing predicted DO to an independent set of DO measurements conducted 1-2 days after initial sampling (10 sites per sample date) using linear regression. We found a strong association between predicted and actual DO values (r 2 = 0.92, n = 120; P \ 0.01), indicating that maps generated in ArcGIS provided good approximations of the spatial variation in DO present in the lake at each sampling occasion. Finally, we estimated the proportional volume of the hypoxic zone for each sampling occasion by calculating the proportion of the upper 1 m of the water column with DO B 4 mg/l.
We used multiple linear regression to assess relationships between environmental variables and lake-wide average DO concentrations. For air temperature, we used maximum daily temperature (MDT) as the covariate, postulating that days with the lowest maximum temperature would produce the most severe winter conditions. Initial analysis revealed little evidence for multicollinearity among variables (Zuur et al., 2009). Model comparison was performed using Akaike's information criterion adjusted for small sample size (AICc; Burnham & Anderson, 2002). We assessed relative model support by considering the variables included in the best-supported model, i.e., the model with the lowest AIC c score (Burnham & Anderson, 2002), and by considering all models \ 7 AIC c units as plausible.

Winter severity
Based on the association between MDT below zero and low DO levels in the lake (see Results), we calculated an index of annual winter severity, and ostensibly the extent of winter hypoxia, by the number of days with MDT \ 0°C during ice cover (November to March) for each year of available air temperature data, spanning 28 years (1987 to 2015) from the Red Rock weather station. We then used linear regression to assess if the number of days with MDT \ 0°C has changed across this time span by assessing if the slope of the trend line was different from zero (Pederson et al., 2010;Isaak et al., 2012).

Results
Ice cover across the lake formed by November 15, 2013 and November 25, 2014. Ice thickness increased to * 0.5 m during both winters. The extent of ice cover was uniform across the lake except around intermittent, small openings (\ 100 m 2 ) at the mouth of Elk Springs, Red Rock, Shambow, and Grayling creeks. Ice breakup began in the first week of April in 2014 and the last week of March in 2015. Ice cover was completely absent by mid-April of both years.
Lake-wide average DO concentration prior to ice formation in the fall was 9.0 mg/l and 10.2 mg/l in 2013 and 2014, respectively, and was uniform throughout the lake regardless of depth and site (Fig. 2). Within 40-50 d after the onset of lake ice, lake-wide average DO declined substantially to minimum levels of 4.3 mg/l (0.2 mg/l per day and 5.9 mg/ l (0.1 mg/l per day) during the winters of 2013-2014 and 2014-2015, respectively. After the sharp initial decline in DO following ice-over to the peak lows at midwinter, under-ice DO stabilized at low levels for 1 to 1.5 months, followed by a rapid increase to pre-ice levels by March 27, 2014 andFebruary 16, 2015. During both years, DO rebound occurred when ice still completely covered the lake, approximately 2-4 weeks prior to ice breakup. During rebound, lake-wide average DO increased from 4.3 to 8.7 mg/l in 28 days in 2014 (0.2 mg/l per day), and from 5.9 to 9.8 mg/l in 14 days in 2015 (0.3 mg/l per day).
Despite a maximum water depth of only 2 m, a significant DO and temperature gradient developed in the water column by midwinter during both years. In contrast to the sharp drop in lake-wide average DO, mean surface DO measured at the ice-water interface remained high during all sampling periods (* 6.5 mg/l, 2013-2014; * 8.0 mg/l, 2014-2015), whereas DO decreased sharply with depth, with average DO at 1 m of 2.6 mg/l (2013-2014) and 5.1 mg/l (2014)(2015), and with DO concentrations often anoxic in the lower hypolimnetic layer [ 1 m deep (Fig. 3). The DO gradient persisted throughout the DO minima period, but ended when under-ice concentrations returned to an unstratified state by March 27, 2014 and February 16, 2015. A strong temperature gradient also developed during ice cover. For example, water temperatures at the ice-water interface ranged from 0.3 to 0.6°C, whereas average lake bottom temperatures ranged from 3.6 to 3.7°C (Fig. 4).
Ice thickness increased over the course of both winters (0.4 cm/day in 2013-2014, and 0.2 cm/day in 2014-2015) but varied between study years. Maximum ice thickness was * 43% greater in 2014 (61.7 cm) than in 2015 (43.1 cm) (Fig. 5) but was generally uniform across the lake at each sampling occasion (variation \ 10 cm). Similarly, maximum snow depth was nearly two-fold greater in 2014 (18.1 cm) than 2015 (9.6 cm). Snow depth, however, was variable across the lake due to wind scouring the ice surface (Fig. 5). Snow was largely absent from the ice surface by March 27, 2014, and by February 16, 2015.
DO interpolation maps showed the expansion and contraction of hypoxic areas as well as differences in DO concentration between the two study winters (Figs. 6, 7). Prior to, and shortly after, complete winter ice cover, the lake had uniformly high DO levels. During midwinter when DO reached minima levels, low DO areas were prevalent in the shallow northwest and southern parts of the lake, particularly in the lowest DO winter of 2013-2014. In contrast, areas near inlet streams remained well-oxygenated throughout the midwinter DO minima period. In 2014, the proportional volume of the hypoxic (\ 4 mg/l) area in the upper 1 m expanded from 7% in early-January to 31% in late February, followed by a decline to 0% by the end of March. In contrast, the proportional volume of hypoxic area of the surface layer was * 1% in 2014-2015, and normoxic DO levels ([ 4 mg/l) were present throughout the lake even in mid-February.
Lake-wide average DO over the two study years showed strong positive association with maximum daily air temperature MDT (r 2 = 0.79, P \ 0.01), and a negative association with snow depth (r 2 = 0.68, P \ 0.01) (Fig. 8). We found no association between  (Table 1). AIC rankings indicated that other models that included snow, ice, and water temperature also influenced average lake-wide DO levels, but to a lesser degree. Plots of MDT and lake-wide average DO indicated that winter minima DO periods during both winters were closely associated with periods when MDT was \ 0°C (Fig. 9). Over the 28-year weather record, there was high variability in the number of days within a given year when MDT was \ 0°C, ranging from 62 to 103 d (mean = 81 days; Fig. 10). Study winters, with relatively low and high DO levels, were above and below, respectively, the long-term average winter severity, with the winter of 2013-2014 having 23% more days when MDT \ 0°C than in 2014-2015 (2014: 88 days; 2015: 68 days). The plot of MDT days per year showed a declining but nonsignificant trend over the 28-year time period ( Fig. 10; slope = -0.236, P = 0.36).

Discussion
Our study demonstrated pronounced spatial and temporal variation in DO availability during winter in a shallow, ice-covered lake. During midwinter, hypoxic DO concentrations were widespread, but some areas under the ice, particularly around stream mouths and immediately below the ice surface, remained well-oxygenated throughout the winter. This spatial variation in DO is not typically observed in studies of winter oxygen depletion and winterkill risk, which have generally utilized single point DO sampling to estimate extent of DO depletion in icecovered lakes (Barica & Mathias, 1979;Babin & Prepas, 1985;Prowse & Stephenson, 1986;Leppi et al., 2016). However, the few studies that have employed more extensive DO sampling such as our study have similarly found evidence for high spatial variability in DO, even in small lakes prone to winterkill (Hasler et al., 2009), suggesting that such variability might be more widespread, especially in lakes with surface inflow (Gangloff, 1996;Guenther & Hubert, 1991). The extent of normoxic conditions in Upper Red Rock Lake in the winter was surprising as near-anoxic conditions had been measured in previous work and the lake has all the attributes common to lakes with high winterkill risk. We measured very low DO concentrations after ice formation for an extended period during both study winters, but most hypoxic conditions were observed at depths [ 1 m deep. Extensive hypoxia in the hypolimnion in lakes has been observed in other studies (Mathias & Barica, 1980;Prowse & Stephenson, 1986;White et al., 2008;Leppi et al., 2016). We also observed hypoxia in the epilimnetic surface layer \ 1 m below the ice surface, but its extent was much more variable than reported in most previous studies of winterkill-prone lakes.
Previous investigations of winterkill lakes reported relatively constant DO depletion rates once ice cover was achieved, allowing for prediction of winterkill based on measurements of lake volume and depletion rate (Barica & Mathias, 1979;Prowse & Stephenson, 1986;Leppi et al., 2016but see Guenther & Hubert 1991 for exception). Many of these studies relied on data from comparatively small lakes with little to no input of surface water from tributaries. In contrast, Upper Red Rock Lake had a distinctive non-linear DO depletion rate and considerable year-to-year variation. If the high DO depletion rate we observed in early winter had continued unchanged, average lake-wide DO concentrations would have reached * 0 mg/l by early-February in 2014, and * 4 mg/l by early-March in 2015. Instead, the DO depletion rate leveled off by midwinter and remained stable until the rapid rebound in DO prior to ice-off during spring.
High variation in location and extent of hypoxia within the lake suggests an unidentified input of DO to the lake during the midwinter oxygen minima, preventing DO from being fully depleted. Plausible mechanisms of DO inputs during midwinter include: (1) input of oxygen-rich stream water, (2) photosynthetic production of oxygen by algae and submergent plants, and (3) oxygen input from springs located on the lake bottom. Our observations of intermittent ice cover on inlet streams and snow-free parts of the lake suggest mechanisms 1 and 2 could be responsible for DO inputs during the midwinter DO minima. Although mechanism 3 is possible, we did not find  Fig. 8 Scatterplots of the relationship between lakewide average dissolved oxygen and maximum daily air temperature (A), lakewide average ice thickness (B), lake-wide average water temperature (C), and lake-wide average snow depth (D) Table 1 DAICc, R 2 and model coefficient estimates with 95% confidence intervals for models of lake-wide average dissolved oxygen Notably, we found a negative association between lake-wide DO and snow depth, supporting this hypothesis. Previous videography surveys identified green aquatic plants persisting and photosynthesizing throughout the winter months under the ice during a year of low snow cover (Gangloff, 1996). We believe that mechanism 1 is the most important driver of DO dynamics, as Gangloff (1996) and our study found high DO at inlet stream mouths, suggesting that streams delivered DO-rich water throughout the winter months. We hypothesize that the lower density and temperature of these inflowing tributaries would have allowed this oxygen-rich water to spread out along the ice-water interface without mixing with bottom layers (Bergmann & Welch, 1985). Guenther and Hubert (1991) also found that winter hypoxia was less pronounced in lakes with surface inflow. Another unexpected result of our study was the lake-wide spring DO rebound prior to ice loss (Baehr & DeGrandpre, 2002;White et al., 2008). The typical pattern of lake re-oxygenation in winterkill lakes is DO rebound concomitant with reconnection with the atmosphere during spring ice breakup (Barica & Mathias, 1979). However, we observed rapid reoxygenation at least 2-4 weeks prior to ice breakup when ice cover was still complete. We hypothesize that increased flows due to valley floor snowmelt along with resumption of photosynthetic activity with   increased light penetration due to dwindling snow depths on the lake ice, influenced the rebound in lakewide DO levels during maximum ice extent in earlyspring. Such rebound has not been documented elsewhere to our knowledge and is an important consideration for estimating duration of winterkill risk.
Our results also suggest that variation in maximum daily air temperature is a potentially useful surrogate for estimating the extent and degree of winter hypoxia. Measurements of winter DO levels in two other winters of quite different severity in Upper Red Rock Lake by Gangloff (1996) provide additional support of this relationship. In the winter 1994-1995, DO was at low levels (\ 0.5 mg/l) at lake sampling sites excluding tributary mouths throughout the winter and the number of MDT days below zero was above average (83 days); in contrast, winter DO was much higher throughout the lake the following winter ([ 6 mg/l) and number of MDT days below zero (n = 64) was well below average. If MDT days below zero is representative of winter severity and winter hypoxia as these data suggest, the weather record data and DO measurements among different years indicate a high degree of interannual variability in the threat of hypoxia. Future climate projections indicate less frequent and less severe hypoxic conditions in shallow, ice-covered lakes with warming temperatures during winter (Fang & Stefan, 2000). However, we found only a slight negative decline in winter severity over a 28-year period. Upper Red Rock Lake occurs in a mountain valley with frequent temperature inversions in winter that result in persistent cold temperatures. We suspect that such inversions may buffer some of the effects of rising regional temperatures (Pederson et al., 2010).
Our findings have important implications for measuring the dynamics of winter hypoxia and estimating its impacts on overwintering fishes. In Upper Red Rock Lake, the inlet tributary mouths likely offer DO refugia key for survival of grayling and other species during winters with otherwise widespread hypoxic conditions (Gangloff, 1996;Davis et al., 2019). Identification of such refugia will help conserve these critical habitats (Weber et al., 2013;Tibbles et al., 2018). The approach we used in our study provides a useful framework for quantifying and mapping the seasonal dynamics of the extent and severity of winter hypoxia, and for identifying critical winter habitats in other systems.
In our study, we show evidence of substantial variation in development of hypoxia both within and among years, emphasizing the need for better understanding of the processes that affect hypoxia over space and time. Better defining the relative role of inlet streams versus in situ photosynthesis on driving midwinter DO levels and re-oxygenation prior to  ice-off would be particularly fruitful areas for further research.