Thermal stratification and meromixis in four dilute temperate zone lakes

Four adjacent lakes (Arco, Budd, Deming, and Josephine) within Itasca State Park in Minnesota, USA, are reported to be meromictic in the scientific literature. However, seasonally persistent chemoclines have never been documented. We collected seasonal profiles of temperature and specific conductance and placed temperature sensor chains in two lakes for ∼1 year to explore whether these lakes remain stratified through seasonal mixing events and what factors contribute to their stability. The results indicate that all lakes are predominantly thermally stratified and are prone to mixing in isothermal periods during spring and fall. Despite brief, semi-annual erosion of thermal stratification, Deming Lake showed no signs of complete mixing from 2006–2009 and 2019–2022 and is likely meromictic. However, the other lakes are not convincingly meromictic. Geochemical data indicate that water in Budd Lake, which contains the most water, is predominantly sourced from precipitation. The water in the other three lakes is of the calcium–magnesium–bicarbonate type, reflecting a source of water that has interacted with the deglaciated landscape. $\mathit{\delta }{}^{\mathrm{18}}{\mathrm{O}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ and $\mathit{\delta }{}^{\mathrm{2}}{\mathrm{H}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ measurements indicate the lakes are supplied by precipitation modified by evaporation. Josephine, Arco, and Deming lakes sit in a valley with likely permeable sediments and may be hydrologically connected through wetlands and recharged with shallow groundwater, as no streams are present. The water residence time in meromictic Deming Lake is short (100 d), yet it maintains a large reservoir of dissolved iron, indicating that shallow groundwater may be an additional source of water and dissolved ions. All four lakes develop subsurface chlorophyll maxima layers during the summer. All lakes also develop subsurface oxygen maxima that may result from oxygen trapping in the spring by rapidly developed summer thermoclines. Documenting the mixing status and general chemistry of these lakes enhances their utility and accessibility for future biogeochemical studies, which is important as lake stratification and anoxia are becoming more prevalent due to changes in climate and land use.

1 Introduction

Four lakes in Itasca State Park, Minnesota, USA – Arco, Budd, Deming, and Josephine (Fig. 1) – have been described as meromictic in the scientific literature since the 1960s (Baker and Brook, 1971; Anderson et al., 1985; Stewart et al., 2009). However, there have never been physical and chemical data from all lakes during all seasons – especially spring and fall when mixing would be expected – to document through the maintenance of a chemocline that they are meromictic. The identification of meromictic lakes is important as they are critical analogues for the understanding of the biogeochemistry of past oxygen-stratified oceans (Degens and Stoffers, 1976). In that regard, the four lakes in Itasca are particularly well suited, if they are meromictic, as they have low amounts of sulfur and abundant iron, resulting in anoxic and “ferruginous” conditions under stratification that are chemically similar to oceans during the Precambrian Era (Swanner et al., 2020).

Figure 1 Map of the four study lakes within Itasca State Park, Minnesota. Basemaps from USGS and ESRI.

Meromictic lakes are characterized by a sharp increase in electrical conductivity at a chemocline separating a cold and/or dense monimolimnion from the epilimnion that persists year-round (Boehrer and Schultze, 2008). Thermal stratification is typically not the main mechanism stabilizing meromixis in lakes of the temperate zone as it is in the tropics due to the erosion of the thermocline in autumn (Boehrer et al., 2017). Meromictic lakes are present on all continents but are thought to be rare (Hall and Northcote, 2012; Stewart et al., 2009). However, seasonally persistent stratification and even meromixis are becoming more common as a result of temperature increases driven by climate change (Ficker et al., 2017). Furthermore, increased monimolimnion water density can also result from increased chemical constituents that enrich the water density of monimolimnion, such as those from runoff, groundwater inputs, remineralization processes associated with organic carbon loading (Hakala, 2004), and deicing salt contamination (Koretsky et al., 2012).

While anoxia and hypoxia are common features of temperate lakes under thermal stratification (Nürnberg, 1995), meromictic lakes have a permanently anoxic monimolimnion. These conditions restrict the habitat of aerobic organisms, including zooplankton and fish; can amplify microbial or geochemical cycling of elements across the redoxcline (Busigny et al., 2016); and can impact the emission of greenhouse gases such as methane (Lambrecht et al., 2020). The investigation of meromictic lakes can help to understand the alterations to global biogeochemical cycles that may result from climate change and land-use changes that are driving increases in hypoxia in lakes (Jenny et al., 2016).

Finally, the sediments of meromictic lakes are valuable as archives of past climatic transitions, vegetation changes, changes in sediment transport, and atmospheric deposition patterns (O'Sullivan, 1983). For instance, meromictic Crawford Lake in Ontario, Canada, was recently promoted as a “golden spike” for the proposed Anthropocene Epoch based on radionuclide patterns recorded in its laminated sediments (McCarthy et al., 2023). Laminated sediments that form in many meromictic lakes are also useful for studying the formation, deposition, and diagenetic transformation of chemically precipitated minerals (Ledesma et al., 2023). Such minerals can potentially be used as proxies for paleoredox conditions, which have application for paleo-marine and paleo-lacustrine studies (Swanner et al., 2020).

The goals of this study are to (1) determine whether Arco, Budd, Deming, and Josephine lakes are meromictic and the factors that drive stratification; (2) investigate the water type, sources, and reasons for meromixis; and (3) describe the unique biological features of these lakes. To achieve these goals, we conducted fieldwork from 2006–2009 and 2019–2022 and integrated data from University of Minnesota students working at the Itasca Biological Station and Laboratories since the 1950s. The description of the lakes, their mixing status, and the factors that control it, as well as key biological characteristics, will provide an expanded understanding of why lakes resist mixing and highlight the utility of these lakes for further scientific investigations.

2 Materials and methods

Formation of the four study lakes occurred during the late Wisconsin glaciation ∼12 000 years ago (Marine Isotope Stage 2; Jennings and Johnson, 2011). They occupy a tunnel valley that was formed beneath the Wadena lobe of the Laurentide Ice Sheet. Following glacial retreat, the melting of stagnant ice blocks within the tunnel valley left depressions in the landscape now occupied by lakes and wetlands (Wright, 1993).

Today, the four lakes investigated in this study sit in the HUC-12 watershed that sources the headwaters of the Mississippi River (Watershed Boundary Dataset, 2022). Budd (478.6 m above mean sea level; m a.m.s.l.) lies at the highest elevation, while Arco (465.8 m a.m.s.l.) and Josephine (465.4 m a.m.s.l.) lie at similar elevations (Fig. 1). Deming has the lowest elevation of the lakes (464.8 m a.m.s.l.).

Lake depth measurements were collected using a Garmin STRIKER 4 dual-beam transducer (sonar) attached to a rowboat or canoe. Depth and GPS measurements were taken every 6 s while the boat was in motion. A Garmin GLO 2 GPS receiver and the ArcGIS Collector app were used to navigate and track the boat's course and ensure even coverage. The shoreline of the lakes was obtained by walking along accessible areas of the shore with the Garmin GLO 2 GPS receiver or from lidar-derived digital elevation models. Bathymetry rasters (1 m resolution) were generated from the depth measurements in ArcGIS Pro 3.0 using a third-degree local polynomial interpolation. These rasters were used to calculate lake volumes and contour maps. Rasters and volume data have been deposited with the Environmental Data Initiative (Swanner et al., 2023).

Morphometric data were used to calculate relative depth (Z_r; %), which is the ratio of the maximum lake depth (Z_m; m) to the diameter of a circle of area equal to that of the lake (A₀; m²), expressed as a percentage in Eq. (1) (Wetzel, 2001):

$$\begin{array}{}\text{(1)}& Z_{\mathrm{r}}={\displaystyle \frac{Z_{\mathrm{m}}}{\mathrm{2}}}\sqrt{{\displaystyle \frac{\mathit{\pi }}{A_{\mathrm{0}}}}}\times \mathrm{100}\end{array}$$

Chemical, physical, and biological parameters measured on the four lakes included depth, temperature, specific conductance, salinity, turbidity, pH, oxidation–reduction potential, dissolved oxygen, photosynthetically active radiation, chlorophyll a, and phycocyanin. Major cations, anions, and isotopes of water ($\mathit{\delta }{}^{\mathrm{2}}{\mathrm{H}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ and $\mathit{\delta }{}^{\mathrm{18}}{\mathrm{O}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$) were determined on lake water retrieved from different depths within the four lakes. Taxon-specific chlorophyll a fluorescence was collected with a FluoroProbe (bbe Moldaenke). The data and description of methods are available from the Environmental Data Initiative (Swanner et al., 2023). Measurements were made and samples were collected from a boat anchored within the deepest basin of each lake.

A string of temperature loggers (HOBO Water Temp Pro v2) placed at different depths was deployed into the deep areas of Deming, Arco, and Budd lakes for 1 year. A conductance logger (HOBO Fresh Water Conductivity Data Logger) was added near the bottom of the strings in Arco and Budd after 6 months. These sensors measured temperature every 30 min and specific conductivity every 2 h. The sensor string was not retrieved from Budd, as it could not be located in May 2022. Conductance measurements with a Yellow Spring Instruments (YSI) ProDSS temperature and conductivity sensor on deployment and removal were used to check for drift in the HOBO conductance logger. Hourly wind speed data for the duration of sensor deployment utilized the ITCM5 (47.2400, −95.1900; 451.1 m elevation) weather station in Itasca State Park. Data were downloaded from MesoWest (https://mesowest.utah.edu/, last access: 12 October 2022). Plots and analyses were produced in Python or R Studio (build 2022.07.2).

The density, Brunt–Väisälä or buoyancy frequency (N², s⁻²), and Schmidt stability were calculated using data from seasonal profiles or the loggers in the RLakeAnalyzer package v.1.11.4.1 with equations described in that package's documentation (Winslow et al., 2019).

Major anions (${{\mathrm{CO}}_{\mathrm{3}}}^{\mathrm{2}-}$, ${{\mathrm{HCO}}_{\mathrm{3}}}^{-}$, Cl⁻, and ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$) and cations (Na⁺, K⁺, Ca²⁺, and Mg²⁺) were used to produce a Piper diagram in Geochemist's Workbench 15.0. The concentration of cations and anions was calculated as the percentage of total cations and anions in milliequivalents (meq) L⁻¹.

The isotopes of water ($\mathit{\delta }{}^{\mathrm{2}}{\mathrm{H}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ and $\mathit{\delta }{}^{\mathrm{18}}{\mathrm{O}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$) were measured on spring or seep water that had been filtered with 0.45 μm nylon syringe filters and stored at 4 °C with minimal headspace until analysis. Samples were analyzed with a Picarro L1102-i isotopic liquid water analyzer at the Stable Isotope Laboratory at Iowa State University. The analytical uncertainty and average correction factor for $\mathit{\delta }{}^{\mathrm{18}}{\mathrm{O}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ are ±0.05 ‰ and ±0.30 ‰ for $\mathit{\delta }{}^{\mathrm{2}}{\mathrm{H}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ relative to the Vienna Standard Mean Ocean Water (V-SMOW).

Samples for microscopy and water color were collected in amber bottles with a Van Dorn sampler from three different depths in each lake, including the subsurface chlorophyll maximum layer (SCML), if present, as determined with the YSI ProDSS. Water color was determined on water filtered through a GF75 glass fiber filter (Advantec) (Cuthbert and del Giorgio, 1992). Absorbance was measured at 440 and 750 nm. The absorption coefficient (g; m⁻¹) was calculated by subtracting the absorbance at 750 nm from the absorbance at 440 nm and dividing the result by the path length (m):

$$\begin{array}{}\text{(2)}& g_{\mathrm{440}}=(\mathrm{2.303}\times A_{\mathrm{440}}-A_{\mathrm{750}})/\left(\text{path length}\right).\end{array}$$

A conversion was necessary to determine the color (mg Pt L⁻¹) of the lake water:

$$\begin{array}{}\text{(3)}& \text{color}=\mathrm{18.216}\times \left(g_{\mathrm{440}}\right)-\mathrm{0.209}.\end{array}$$

Water sampled from the SCML was preserved with 1 % Lugol's solution upon returning to the laboratory. Fixed samples were settled in the dark for 3 to 7 d.

Student reports from courses taking place over several decades at the Itasca Biological Station and Laboratories (IBSL) (Knoll and Cotner, 2018), formerly the Itasca Biological Station, were acquired from the library at the University of Minnesota, Twin Cities. Morphometric data, temperature and conductance profiles, observations of water column mixing, and inferences from biological experiments in these reports are referenced in the current study.

3 Results and discussion

To ascertain if the study lakes are meromictic, morphometric data were acquired to assess whether the basin shape was prone to stratification (Sect. 3.1). Specific conductance profiles collected in this study and from historical records are used to determine if seasonal mixing occurred. Where those records are lacking or incomplete, temperature profiles or high-frequency sensor data are used to assess the stability of stratification. Because of the importance of specific conductance measurements in indicating mixing, the major ions that contribute to the lakes' conductance were determined (Sect. 3.2). These data along with watershed characteristics are also used to assess potential water and ion sources to the lakes. Finally, because meromixis can result in distinctive biological zones and biogeochemical reactions therein, an assessment of these major features was made and compared to historical biological records (Sect. 3.3).

3.1 Morphometry and mixing status

All four study lakes exceed a Z_r of 4 % and have surface areas of less than 500 000 m², typical of meromictic lakes in the temperate zone (Swanner et al., 2020). Deming Lake is the largest and deepest of the four study lakes with a 54 325 m² surface area and a maximum depth of 20.8 m (Table 1), with 4.3 % of the surface area encompassing water depths at or below 17 m. The maximum depth of Deming Lake has previously been reported as 16.5 m (Hooper, 1951), 17 m (Baker and Brook, 1971), and 17.6 m (Lascu and Plank, 2013). Josephine Lake has the longest fetch and lowest Z_r of all study lakes. The maximum depth of Josephine Lake (14.8 m) has previously been reported as 10.3 m (Baker and Brook, 1971), 12–13 m (Callis et al., 1976), and 13 m (Gage and Gorham, 1985). The surface area with depths at or below 13 m in Josephine Lake is 5.1 %. Budd Lake has the highest Z_r (8.8 %). The greatest variation was between the previously reported maximum depth of Budd Lake (10.8 m; Baker and Brook, 1971) and measurements presented here (16.1 m). Arco Lake has the smallest surface area of 24 180 m² and a maximum depth of 12.8 m. Arco and Budd lakes have steep banks along their eastern sides (Fig. 1).

Table 1 Morphometric parameters for the four study lakes.

In all lakes, steep-sided deep holes were detected that were deeper than previously reported (Figs. 2 and S1 in the Supplement). Differences in technology and the small surface areas encompassing the deepest parts of the basins could account for the variation – these basins may have been difficult to locate and map with manual techniques. Alternately, the temperature, pressure, and/or salinity gradients present at the tops of these holes could deflect sound beams, causing errors in our depth measurements (Boehrer and Schultze, 2008). Maximum depths should be verified with measurements by multiple techniques in the future.

Figure 2 Bathymetric maps of the study lakes overlain on the digital elevation models. Depth is in meters. The fence diagrams correspond to the cross sections in Fig. S1.

Deming Lake has the longest record of seasonal profiles of temperature and specific conductance (2006–2009 and 2019–2022), as its sediments have previously been used to investigate Holocene climate variations (Lascu et al., 2012; McLauchlan et al., 2013). A chemocline, or sharp increase in specific conductance, persists in Deming Lake at all sampled times (Fig. 3), which is expected for a meromictic lake. The chemocline occurs at 11–13 m (Fig. 3). The chemocline was below 10 m in 1989 (Church et al., 1989) and from 11–13 m a decade later (Reiter et al., 1998) based on temperature-compensated conductance measurements from IBSL student reports, indicating that the chemocline has been at a consistent depth for decades.

Figure 3 Specific conductance profiles of (a) Arco Lake, (b) Budd Lake, (c) Josephine Lake, and (d) Deming Lake. For Deming Lake, seasons were classified according to solstice and equinox dates. Note the variability in the scales of x and y axes.

A systematic increase in the magnitude of specific conductance readings below the chemocline in Deming Lake was observed between the 2006–2009 and the 2019–2022 datasets (Fig. 3). This could result from variations in the amount of in-lake primary productivity over time (Campbell, 1977) or could increase with time due to a mixing event (Katsev et al., 2010). Deming Lake mixed after a beaver dam broke on the western side of the lake in 1997 (Frane and Walberg, 1997). If the lake has not mixed since aside from this catastrophic event, the specific conductance should have increased from 1997 to 2006–2009 and from 2006–2009 to 2019–2022 at similar rates. The 1997 conductivity data were converted to specific conductance using the reported temperature values (Frane and Walberg, 1997), and the value at 14 m in July 1997 was compared to the values at 14 m in July 2006 and July 2021. The specific conductance values increased during both intervals, but the rate changed from 31 $\mathrm{\mu }\mathrm{S}{\mathrm{cm}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ from July 1997 to July 2006 to 11 $\mathrm{\mu }\mathrm{S}{\mathrm{cm}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ from July 2006 to July 2021. This could imply a partial mixing event in the years between 2006 and 2021 or some change to the supply of water and dissolved ions to the lake between the two time intervals being compared.

Changes in precipitation could have affected the rate of the increase in specific conductance in Deming Lake over this period. Historical drought records are available for the region encompassing Itasca State Park from the year 2000 (National Integrated Drought Information System, 2022; Fig. S2 in the Supplement). Drought over this period could have modulated the rate of the lake's specific conductance increase through increased evaporation (Jellison and Melack, 1993). Alternatively, or in addition, the magnitude of water sources of differing specific conductance supplying the lake may have been altered (Ludlam and Duval, 2001). In the second case, drought would decrease the magnitude of low specific conductance precipitation and increase the proportion of higher specific conductance groundwater to the annual water budget. As Deming Lake does not have a permanent surface inlet, changes in streamflow can be neglected. The amount of recharge to the groundwater system may have also decreased due to drought, effectively decreasing groundwater inputs to the lake and counteracting the effects of evaporation on specific conductance in Deming Lake.

The record of seasonal specific conductance profiles, which is necessary to document the seasonally persistent chemocline that characterizes meromictic lakes, is only available for the period of 2019–2022 in the other three lakes. We therefore take a comparative approach to establish the likely mixing status of the other lakes. Budd Lake had the lowest range of specific conductivity values (23–90 µS cm⁻¹). A weakly demarcated chemocline was observed around 4–5.5 m water depths for May 2021, July 2021, October 2021, and May 2022. The January 2021 profile shows a deeper chemocline (Fig. 3). Specific conductance values were uniform down to 8 m in Josephine Lake during the current study period (Fig. 3). A chemocline was only detected below 9 m in October 2021, while depths below 11 m were not assessed on other dates during the study period due to the difficulty of finding the deepest area before mapping in 2022. Autumnal circulation down to 10 m was reported in Josephine Lake in 1975 (Gage and Gorham, 1985). A chemocline was also present at all time points for Arco Lake between 7–10 m except for May 2022 (Fig. 3). It was poorly developed in May 2021, and its deepest occurrence was in January 2021. Arco Lake has a maximum depth of 12.6 m, but no specific conductance data were collected below 12 m as part of the profiles acquired during the study period due to the difficulty of finding the small deep spot (representing 1.7 % of the surface area). The chemocline is shallowest in mid-summer, which is consistent with the 7 m chemocline observed in temperature-equilibrated conductivity readings in 1975 and 1976 (Evans and Bjerklie, 1975; Barnes et al., 1976) and 2011 (Harren et al., 2011).

Arco Lake has the next best documentation of chemocline dynamics after Deming Lake due to the placement of temperature sensors (1.5, 3.5, 5.5, 7.5, 9.1, 9.5 m) and a conductance sensor (9.1 m) from May 2021 to May 2022 (Fig. 4). In the summer of 2021, Arco Lake developed a summer thermocline. The temperatures at different depths started to converge in late fall until the lake became isothermal and continued cooling from 6 to 4 °C in mid-November 2021, representing the erosion of thermal stratification. The Brunt–Väisälä or buoyancy frequency (N), commonly expressed as N² in units of s⁻², calculated from the time series temperature data is near zero at that time, indicating an erosion of stability, which persisted until the sensors were removed on 16 May 2022, over a week after the ice-off on 8 May (Fig. S3 in the Supplement). A drop in the conductance value at 9.1 m in early November 2021 preceded isothermal conditions and indicated mixing down to at least 9.1 m (Fig. 4). The dimensionless lake number was calculated according to Imberger and Patterson (1989) using the time series and hourly wind data from Itasca State Park (Fig. S4 in the Supplement). It started to consistently rise on 15 May 2022, as the epilimnion warmed. In combination with the seasonally variable chemocline depth (Fig. 3), these analyses indicate that Arco Lake mixes during isothermal periods and is not a meromictic lake. The data presented are most consistent with Arco being a holomictic lake.

Figure 4 (a) HOBO sensor data for temperature and conductance from Arco Lake, collected from May 2021 to May 2022 and from October 2021 to May 2022, respectively. (b) Isotherm plot of the temperature data from the HOBO sensor.

Data from the temperature sensors (0.5, 2.5, 5.5, 8.5, 11.5, and 14.5 m) in Deming Lake from June 2019 to May 2020 are presented in Fig. S5 in the Supplement. The lake became isothermal in early November 2019. The Brunt–Väisälä frequency calculated from this data dropped to near zero and rebounded slightly due to under-ice thermal stratification, a phenomenon not observed in Arco Lake (Fig. S3). The specific conductance profiles for Deming consistently show the same trend, although the exact depth and magnitude of the chemocline vary seasonally and on the decadal scale (Fig. 3). This indicates partial seasonal mixing in Deming, likely during isothermal periods, but it may be insufficient or of too short a duration to fully mix the lake. Deming was ice-free on 25 April 2020 (Lake Ice Out Dates, 2022), and the lake number started rising immediately due to the onset of thermal stratification (Fig. S4), limiting the opportunity for mixing. Long-term seasonal observation of a chemocline is consistent with Deming being a meromictic lake (Zadereev et al., 2017).

Arco therefore presents an end-member of a holomictic lake and Deming an end-member of a meromictic lake among the four lakes studied here. In the absence of high-frequency sensor profiles for Budd and Josephine lakes, the Schmidt stability (Fig. S6 in the Supplement) and the Brunt–Väisälä frequency (Figs. S7–S10 in the Supplement) were calculated from seasonal profiles of temperature (Fig. S11 in the Supplement) and salinity for all four lakes. The similar volumes, depths, and surface areas of the lakes (Table 1) permit comparison of the Schmidt stability results. The Schmidt stability indicates that Deming and Budd lakes are the most strongly stratified, followed by Josephine and Arco lakes. It is unknown when ice came off Budd Lake in 2022, but Deming Lake was ice-free on 7 May, and Arco and Josephine lakes were ice-free on 8 May (Lake Ice Out Dates, 2022) following strong winds on the evening of 7 May. Likely, Budd Lake was also ice-free by 8 May given its similar morphometry and location. A chemocline consistent with other summer observations was present in Budd Lake on 17 May 2022, which could indicate a lack of spring mixing. However, the January 2021 specific conductance profile for Budd Lake shows a very weak and deep chemocline, which would be expected following an autumn mixing event. Without full profiles of specific conductance in Josephine Lake, it is impossible to evaluate if the lake mixed. However, based on existing profiles, if a chemocline persists through spring turnover in Josephine Lake as suggested by the sharp and deep chemocline observed in May 2022, then the monimolimnion must be limited to depths below 11 m, representing only 3.6 % of the total lake volume (11 404.5 m³ of 313 295.5 m³).

Water color is related to the abundance of dissolved organic carbon (DOC) (Pace and Cole, 2002). Of the four study lakes, Budd Lake had the most colored water, followed by Deming Lake (Table S1 in the Supplement). Arco and Josephine lakes had very little color. Visually, Budd and Deming lakes appear brown, while Arco and Josephine lakes appear green. Enhanced water color could lead to a shallower thermocline and stronger stratification in Budd and Deming lakes due to greater light absorption by compounds conferring color (Houser, 2006). The Brunt–Väisälä frequencies, calculated from temperature and salinity profiles, were generally highest in the epilimnia during late summer (July 2021; Figs. S7–S10). This indicates that temperature is more important than salinity to the stability of stratification in these lakes. In August 2022, when Deming Lake had the highest N² values observed from the 2019–2022 dataset (data not plotted in Fig. S9), the meromictic stability (S^′) calculated as described by Walker (1974) was 7.81 J m⁻², within the range observed for the meromictic and ferruginous Lake 120 in Canada (Campbell, 1977) and Lake Nordbytjernet in Norway (Hongve, 1999). These observations suggest that meromixis in Deming Lake is primarily due to thermal stratification, with a small contribution from salinity that prevents mixing in this lake but maybe not in the others. This phenomenon could be common to dilute meromictic lakes. The mixing event at Deming Lake in the summer of 1997 due to a breached beaver dam on the western side of the lake (Frane and Walberg, 1997) indicates that mixing can occur due to catastrophic events.

3.2 Chemical characteristics and water sources

A Piper diagram was used to define the water types of the four lakes based on their major cations and anions. The predominant cation in all four study lakes was calcium (Fig. 5). Budd Lake had a greater proportion of milliequivalents from sodium and potassium than the other three lakes, with the charge balance being mostly from magnesium. Bicarbonate and carbonate ions represented nearly all anion milliequivalents in Deming, Arco, and Josephine lakes, whereas chloride and sulfate contributed to the anion balance in Budd Lake. The water type for Deming, Arco, and Josephine lakes is a calcium bicarbonate type, whereas Budd Lake does not have a dominant water type. The low specific conductance and lack of water type in Budd Lake indicate the water source is predominantly precipitation. Calcium, magnesium, and carbonate–bicarbonate ions are sourced from the calcareous Itasca Moraine and are typical of the Itasca region (Megard et al., 1993). This indicates that water supplying Deming, Arco, and Josephine lakes has undergone more water–rock reaction than the water supplying Budd Lake, such as might be expected from an increasing proportion of groundwater seepage to a lake's water inputs. Deming, Arco, and Josephine lakes lie in a tunnel valley channel containing coarser sands and gravel, while Budd Lake is poorly integrated hydrologically with the other three lakes, and the net seepage in Budd Lake is probably the lowest among all the lakes. Arco Lake (465.8 m a.m.s.l.) is located 51 m north of Josephine Lake (465.4 m a.m.s.l.) and Deming Lake (464.8 m a.m.s.l.) is located 284 m north of Arco Lake. The similar lake levels and the occurrence of wetlands between Arco and Deming lakes (Figs. 1 and S1) indicate these lakes are likely to be hydrologically connected. If the water table mimics the surface topography, groundwater potentially flows from Arco–Josephine to Deming Lake. However, Budd Lake (478.6 m a.m.s.l.) is perched highest in the watershed, located off-axis of Josephine–Arco–Deming lakes, to the west of a ridge (Fig. 1). Because Budd Lake is perched, the ridge is likely composed of less permeable material, possibly till, than the valley sediments. The presence of wetlands on the north end of Budd Lake that follow a small valley towards Deming Lake suggests that there could be a periodic surface or near-surface hydrological connection from Budd Lake to Deming Lake. Lakes highest in their watershed are more likely to receive a greater proportion of water from precipitation (Kratz et al., 1997), which is consistent with the dilute nature of Budd Lake water.

Figure 5 Piper diagram of major cations and anions collected in May 2021.

The isotopic composition of water (i.e., $\mathit{\delta }{}^{\mathrm{2}}{\mathrm{H}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ and $\mathit{\delta }{}^{\mathrm{18}}{\mathrm{O}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ in ‰) in May 2021 varied little with depth in each lake (Fig. S12 in the Supplement). In Arco, Budd, and Deming lakes sampled in June 2019, the epilimnion waters were depleted relative to the deeper water. Additional measurements at Deming Lake in July and October 2021 show enriched values in the epilimnion relative to the deeper water. Preferential evaporation of light isotopes during the summer and early fall drives the epilimnia to more enriched values. In Fig. 6, all lakes' data points diverge from a local meteoric water line (LMWL; Stelling et al., 2021), falling on a lake evaporation line (LEL) with the equation ${\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\left(\mathrm{‰}\right)=\mathrm{4.5}\times {\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\left(\mathrm{‰}\right)-\mathrm{29.9}$. The intersection of the LMWL and the LEL is the composition of regional groundwater, representing isotopically depleted snowmelt as the predominant source of recharge to the local aquifers (Krabbenhoft et al., 1994).

Figure 6 Isotope data (points) and fit (dashed line) of an evaporation line ($\mathit{\delta }{}^{\mathrm{2}}{\mathrm{H}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}=\mathrm{4.5}\times \mathit{\delta }{}^{\mathrm{18}}{\mathrm{O}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}-\mathrm{29.9}$, r-squared=0.93) from the four study lakes plotted with a local meteoric water line (solid; Stelling et al., 2021). Crosses are spring or bog water data from Itasca State Park (Table S2). The cross closest to the lake data points is the Deming bog sample, and the cross closest to the intersection of the lake evaporation line (LEL) and the local meteoric water line (LMWL) is the Nicolet Creek sample. Analytical precision is within the symbol sizes.

To ascertain the isotopic composition of groundwater in the area, historical $\mathit{\delta }{}^{\mathrm{2}}{\mathrm{H}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ and $\mathit{\delta }{}^{\mathrm{18}}{\mathrm{O}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ data from springs at Elk Lake in Itasca State Park were retrieved from the Minnesota Spring Inventory (Minnesota Department of Natural Resources, 2022). The Elk Lake springs were sampled again in May 2022, and samples were collected from spring along Nicolet Creek and the bog on the southeast side of Deming Lake (Table S2 in the Supplement). These sites had abundant marsh marigolds (Fig. S13 in the Supplement), which bloom in May at groundwater discharge sites in northern Minnesota (Rosenberry et al., 2000). The sites at Elk Lake and Nicolet Creek had visible iron mineralization, indicating reducing, iron-bearing water discharge at the surface (Fig. S13). The Nicolet Creek spring sample lies closest to the intersection of the LEL and LMWL, whereas the Deming bog sample lies closest to the lakes but is more enriched than the lakes (Fig. 6).

The seasonal amplitude of a lake's oxygen isotopes can be used to calculate a lake's water residence time (Engel and Magner, 2019). A time series of temperature data from the weather station ITCM5 was used to determine average daily temperatures, and seasonal maximum and minimum temperature values were used to calculate the amplitude in $\mathit{\delta }{}^{\mathrm{18}}{\mathrm{O}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ inputs of precipitation and water vapor $\mathit{\delta }{}^{\mathrm{18}}{\mathrm{O}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ (Fig. S14 in the Supplement). These two values were used to estimate a mean water residence time of 100 d for Deming Lake. The Deming Lake $\mathit{\delta }{}^{\mathrm{18}}{\mathrm{O}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ values encompass early-spring, mid-summer, and late-summer data, but data from the other three lakes were only from May and June time points (Fig. S12), and so the mean water residence time could not be calculated by this method.

Potential water sources to Deming Lake are precipitation, outflow from the boggy areas to the south and west, and groundwater. In addition to higher specific conductance values than the other lakes, Deming also has higher concentrations of iron in the monimolimnion than the other lakes, exceeding 1 mM, and has sediments enriched in biogenic magnetite (Lascu et al., 2012). As springs in Itasca State Park have visible iron mineralization (Fig. S13), it is possible that groundwater supplies iron into Deming Lake (Swanner et al., 2020). Due to active iron redox cycling, which scavenges iron into the monimolimnion in meromictic lakes, dissolved iron can accumulate to high concentrations in its ferrous form (Busigny et al., 2016; Campbell, 1977). The low abundance of sulfur (Fig. 5) means very little sulfide can be produced in the monimolimnion. The lack of anion for precipitation of ferrous iron as sulfide mineral allows for the maintenance of these high ferrous iron concentrations and ferruginous conditions.

3.3 Biological characteristics

An SCML often develops within the metalimnion, the middle seasonally mixed layer of meromictic lakes (Baker and Brook, 1971; Camacho, 2006). Stable stratification, the availability of both light and nutrients, and avoidance of zooplankton grazing in the epilimnion drive the formation and location of the SCML in meromictic lakes (Klausmeier and Litchman, 2001).

Prior work on these four lakes indicated the presence of a lake-wide turbidity maximum layer that was persistently present just below 5 m in Deming Lake (Baker and Brook, 1971). This peak was observed between late May and August from 1968–1970 but was absent in the winter. It was reported in that study that bacteria and phytoplankton both contributed to the turbidity maximum based on microscopic observations. This work indicated the layer was populated by the Cyanobacteria Oscillatoria agardhii var. isothrix (this genus is now called Planktothrix).

From 2019 to 2022, chlorophyll a measured via multi-wavelength fluorescence acquired with a FluoroProbe allowed for the deconvolution of the fluorescence signal so chlorophyll a could be attributed to one of four taxonomic groups: Cyanobacteria, Chlorophyta, Diatoms and Dinophyta, or Cryptophyta. Taxon-specific concentrations of chlorophyll a in each lake inform phytoplankton community structures with depth (Fig. 7). Deming Lake has a persistent SCML around 5–6 m attributable to Cyanobacteria. The SCML was less distinct in May 2022, 1.5 weeks after ice-off. It was also less distinct in May 2021. This suggests that the SCML in Deming Lake is a summer phenomenon that develops at or below the thermocline (Fig. S11), which is consistent with the observations from Baker and Brook (1971). In June 2019, a profile of photosynthetically active radiation was also recorded (Fig. S15 in the Supplement) and showed an inflection point at 5.5 m, corresponding to the depth of the SCML at that time. These observations are consistent with Cyanobacteria forming a dense accumulation at that depth, attenuating the light flux into deeper waters.

Figure 7 Multi-wavelength chlorophyll fluorescence (i.e., FluoroProbe) measurements at (a) Arco, (b) Budd, (c) Josephine, and (d) Deming lakes, showing seasonal and depth trends in major taxonomic groups. Note that the x scales vary between plots.

Samples from the epilimnia or the SCML were taken from the four lakes in May 2022: Arco (5 m), Budd (4.5 m), Deming (3.5 m), and Josephine (4.5 m). The most common phytoplankton of the four lakes were the Cyanobacteria Planktothrix sp. and Aphanocapsa sp. and the green alga Monoraphidium sp. (Fig. S16 in the Supplement). Some other phytoplankton species identified were the red alga Cryptomonas sp., the Cyanobacteria Planktolyngbya sp. and Dolichospermum sp., and other colonial green algae and Cyanobacteria species.

All four lakes have subsurface oxygen maxima which exceed air saturation during the summer but lack this feature in the fall (Fig. 8). The subsurface oxygen maximum in Deming Lake usually occurs at 4 m, near the top of the thermocline (Bieter et al., 1991; Church et al., 1989; Balk et al., 2007). Metalimnetic oxygen maxima can be a result because dissolved oxygen is more soluble in cold water after spring mixing, but as the water warms, the gas solubility decreases, causing oxygen to become supersaturated. Stratification induced by the thermocline then limits the discharge of supersaturated oxygen to the atmosphere (Wilkinson et al., 2015). Deming Lake rapidly develops a thermocline after ice-off (Fig. S5), providing a mechanism for gas trapping. Enhanced biological productivity can also contribute to the oxygen maximum but is generally a smaller contributor (Wilkinson et al., 2015; Craig et al., 1992). The subsurface oxygen maximum is above the SCML in Deming Lake, suggesting that the entrainment of spring oxygen is a larger contributor than photoautotrophy to the oxygen maximum. However, there is also a pH maximum at 4 m during the summer (Bieter et al., 1991; Fig. S18 in the Supplement; Church et al., 1989; Barkow and Habedank, 1990), as is expected for oxygenic photosynthesis. Itasca Biological Station and Laboratories student experiments that quantified photosynthesis and respiration via dissolved oxygen measurements in bottle experiments generally found that net oxygen release was minimal in the SCML due to vigorous respiration (Barkow and Habedank, 1990; Engstrom et al., 1974). However, primary productivity is often highest within the SCML of other meromictic lakes (Camacho, 2006). Measurements of the O₂ : Ar ratio would help to quantify the contribution of oxygen trapping vs. oxygenic photosynthesis to the Deming Lake subsurface oxygen maximum (Craig et al., 1992).

Figure 8 Dissolved oxygen measurements at (a) Arco, (b) Budd, (c) Josephine, and (d) Deming lakes, showing seasonal and depth trends.

The occurrence of a summer SCML in Arco, Budd, and Josephine lakes was more variable than in Deming Lake (Fig. 7). Arco Lake had a spring bloom of Cyanobacteria in May 2021 and in May 2022, at which time the highest presence of chlorophyll a of this study was recorded. At both times, the SCML occurred below the oxygen maximum in a hypoxic zone of the water column. The May 2022 SCML at 5 m in Arco Lake also attenuated light below that depth (Fig. S15), as was observed in Deming Lake in June 2019.

Arco, Budd, and Josephine lakes also have subsurface oxygen maxima during summer, although not as consistently as Deming. When the depth of the thermocline is plotted against the depth of the oxygen maximum for all four lakes, there is a significant correlation if October 2021 samples are excluded (Fig. S19 in the Supplement). In the fall, the subsurface oxygen maximum is absent (the highest dissolved oxygen values are at the surface), and the thermocline is deepest (Fig. 8). This may reflect that in fall the thermocline has deepened into hypoxic or anoxic waters.

4 Conclusions

Although reported to be meromictic, based on seasonal temperature and specific conductance profiles, Arco, Budd, and Josephine lakes may be holomictic. Higher-frequency specific conductance profiles in multiple seasons during turnover events would better establish their frequency of mixing. Deming Lake has been meromictic in the periods of data collection presented here (2006–2009 and 2019–2022) but experienced a mixing event in 1997 when an adjacent beaver dam broke. Although a seasonally persistent chemocline occurs in Deming Lake, stratification is predominantly a thermal phenomenon, with stability mostly conferred by a thermocline that develops rapidly after isothermal spring and fall periods. Thermal stratification is not typical of meromictic lakes in the temperate zone, but this weak stratification may be common in dilute boreal lakes (Meriläinen, 1970; Campbell, 1977) and may be enhanced by water color. However, such meromixis may be easily perturbed by hydrographic or land-use changes (Hongve, 2002), which is consistent with Deming Lake mixing in 1997.

The chemical composition of the study lakes reflects their position in the watershed and water sources. The highest-elevation lake, Budd Lake, was very dilute and had no distinguishable water chemical type and likely received most of its water from precipitation. Arco, Deming, and Josephine lakes are all calcium–magnesium–bicarbonate-type waters (Fig. 5). Deming Lake had the highest ionic load, which is consistent with a greater contribution of water that had undergone some water–rock interaction, likely groundwater. These results further support the conclusion that due to the dilute nature of the water in these lakes, stratification is predominantly a thermal rather than chemical phenomenon. Further work is necessary to characterize groundwater composition and quantify and localize groundwater input to determine if Deming Lake could be classified as crenogenic meromictic. The short water residence time in Deming Lake implied by the analysis of $\mathit{\delta }{}^{\mathrm{2}}{\mathrm{H}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ and $\mathit{\delta }{}^{\mathrm{18}}{\mathrm{O}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ along with the minimal surface water inputs and likely high permeability of tunnel valley deposits suggest the shallow groundwater system should be considered a major water source.

All four lakes have an SCML during summer, as was reported over decades of observations, which is consistent with all four lakes being thermally stratified in the summer months. The depth of the SCML varies in Arco, Budd, and Josephine and is sometimes absent but persists at around 5 m in Deming, corresponding to the base of the photic zone. Subsurface oxygen maxima which exceeds air saturation likely develop by trapping gas in cold water that then warms with rapid thermocline development in the spring. Corresponding pH maxima in Arco and Deming indicate that oxygenic photosynthesis could also contribute to the subsurface oxygen maxima and indicate that zones of maximum oxygen production, primary productivity, and biomass may be spatially decoupled in these lakes due to the strong gradients resulting from stratification.