What killed Frame Lake? A precautionary tale for urban planners

Study area

Frame Lake is located adjacent to the downtown area of the City of Yellowknife (Fig. 1). The lake is eutrophic, and has arguably become endorheic in recent decades, due to disruption of inflow and outflow streams as a result of urban land use change. The lake has a maximum depth of 4.5 m in the Southern Basin and 5.5 m in the northern basin. The Yellowknife area experiences about 100 frost-free days a year, with ice normally covering Frame Lake by the end of October through late May or early June. The maximum depth of ice-cover in small Yellowknife area lakes, including Frame Lake, is approximately 1.0 m (Dirszowsky et al., 2013). The lake is dominated by extensive macrophyte cover, reflecting its highly nutrified state (Fig. S2). The landscape of the Yellowknife area can be described as a low relief (10–20 m) rolling plain, underlain by Canadian Shield geology. Drainage is influenced by this underlying bedrock, with most lakes draining into Great Slave Lake (Dirszowsky et al., 2013).

Anthropogenic impact

The development that would eventually result in the founding of the City of Yellowknife is linked to the discovery of gold in the area in 1896 (Moir et al., 2006; Government of Canada, 2013b). The opening of three major gold mines (Giant Mine, Con Mine and Discovery Mine) in the 1930s led to the rapid expansion and urbanization of Yellowknife through the subsequent decades, with major growth occurring after the city was established as the territorial capital in 1967 (Fig. 2) (Moir et al., 2006). The world class Giant Mine went on to be the most successful mine in the area, producing 7.6 million ounces of gold from the time of first production in 1948 to its closure in 2004. The separation of ore from gold at the Giant Mine was also responsible for the production of 237,000 tonnes of arsenic trioxide (As₂O₃) waste as a by-product of the ore-roasting process (GNWT-RWED, 1993). In the early days of the Giant Mine, a lax regulatory environment resulted in the direct release of As₂O₃ dust into the atmosphere via stack emissions. Although the Giant Mine operators eventually implemented more sophisticated technologies to reduce arsenic aerial emissions by the 1950s, methods which improved over time, approximately 19 million kg of As₂O₃ was released into the atmosphere. These aerial emissions fell on the surrounding landscape, which served as a major contaminant legacy in the area (Palmer et al., 2015; Galloway et al., 2010; Galloway et al., 2017; Government of Canada, 2013b).

The City of Yellowknife eventually grew to encompass the area surrounding Frame Lake. Due to its centralized location within the city and low inflow/outflow rates, the lake is susceptible to sediment loading. Water quality in the lake has thus been hypothesized to have not only been negatively impacted by (1) As₂O₃ originating from aerial mining emissions, but also by (2) the construction of a causeway in the mid-1970s that disrupted outflow, and (3) the construction of roads and built infrastructure in the 1960s and 1970s that disrupted inflow to the lake, and (4) surface runoff of degraded water quality from urbanization around the lake (Fig. 2). There are anecdotal reports of the physical dumping of mine waste and sewage into the lake, and there is documentation that the lake was used as a snow dump, which would have contained a considerable sediment content that would have been scraped up by the plows (D Lovell, pers. comm., 2015). These sediment sources would have contributed to both an increase in sedimentation rate in the lake and would have contributed to the deterioration of the lake’s water quality. An increase in extensive lake bottom vegetation, which upon dying down each winter, has also contributed to an increase in sedimentation in the lake (Dirszowsky et al., 2013; Patterson et al., 2015a). This increase in sedimentation, coupled with low turnover rates and contamination, eventually led to the significant nutrification of Frame Lake and development of winter dysoxic conditions (Dirszowsky et al., 2013; Patterson et al., 2015a).

Field work and sampling design

The field work component of this study was carried out under the Aurora Research Institute licence nos. 14435 and 14965. Three sediment cores were retrieved from the southern basin of Frame Lake using a two-faced rectangular freeze corer and a single faced freeze corer in September 2012 (Fig. 3). The two cores obtained using the two-faced freeze corer were from a depth of ∼4.6 m and were 60.0 cm (core 2FRF1) and 60.4 cm in length (core 2FRF2). The core obtained using the singe-faced corer (1FR) from a water depth of 4.1 m was 86 cm long (see Galloway et al., 2010 for detailed description of core collection method). The freeze cores were shipped to Carleton University for further analysis. The stratigraphy of the cores was logged at Carleton University, and color variation was documented using a Munsell colour chart (Munsell Color Company, Inc, 2000).

Laboratory analyses and chronology

The freeze cores used in this research were subsampled for the various proxy analyses using a custom designed sledge microtome (Macumber et al., 2011) at sampling intervals ranging from 0.1 to 0.5 cm. The single faced core (1FR) was sliced at a resolution varying between 0.2 and 0.5 cm to a depth of 60 cm and the samples were commercially analyzed by ACME Analytical Laboratories, Vancouver, using inductively coupled mass spectrometry (ICP-MS) following Aqua Regia digestion (ICP-MS 1F/AQ250 package) to detect environmentally important metals that are bioavailable (EPA, 2011) and because complete digestion methods, for example multi-acid digestion, can volatilize As, a contaminant of potential concern in this study (Table S1) (Parsons et al., 2012).

Core 2FRF1 was sliced to 0.5 cm resolution to a depth of 15 cm and submitted for radiocarbon dating to the André E Lalonde (AEL) AMS laboratories (AMS, 2014). Additional × range finder radiocarbon dates were obtained from lower in the core. Chronology for the 20th and 21st century uppermost part of the core was based on ²¹⁰Pb dates obtained from a fourth Frame Lake core collected from a few meters away in the same basin in 2014, as well as from three ash layers derived from New Year’s Christmas tree bonfires that were held on the lake in 1968, 1969 and 1971 (V Sterenberg, former Yellowknife Mayor, D Lovell, pers. comm., 2015) (Fig. 3).

Core 2FRF2 was sliced to 0.1 cm resolution for various proxy analyses to a depth of 30 cm. Arcellinida analysis was carried out on subsamples at 1 cm resolution. Three cc’s of wet material from the Arcellinidan aliquot were sieved first at 297 µm then at 37 µm to remove both coarse and fine organic and mineral debris. The samples were then split into six aliquots for micropaleontological analysis using a wet splitter (Scott & Hermelin, 1993). Arcellinida in each sample were quantified using an Olympus SZH dissecting binocular microscope (Olympus, Center Valley, PA, USA) at magnifications ranging between 20 and 50× until counts of at least 150 were obtained. Total carbonates and organic carbon were measured at 1 cm intervals in the 0.1 cm interval immediately adjacent to those used for micropaleontological analysis using loss on ignition (LOI) analysis (Ball, 1964).

Aerial photography

Aerial photos documenting the urbanization around Frame Lake between 1937 and 1996 were acquired through the National Air Photo Library (NAPL)—Natural Resources Canada. Air photos from 1937, 1964, 1976, 1985 and 1996 were used to show major development around the lake (Fig. 2) Google Earth was used to acquire a 2017 image of Frame Lake (Fig. 2).

Statistical analysis

Thirty-two Arcellinidan species and strains were identified using several well-illustrated papers that utilized the strain taxonomic concept as identification guides (e.g., Reinhardt et al., 1998; Patterson et al., 2013; Patterson et al., 2015b; Nasser et al., 2016) in 30 samples and converted to relative abundances (Table 1). For illustrated specimens of the arcellinidan pertaining to local NWT Arcellinida, refer to Nasser et al. (2016). Probable error (pe) was calculated for each sample using the formula: ${\displaystyle pe=1.96\left(\frac{s}{\sqrt{X_i}}\right)}$ where s is the standard deviation of the population and X_i is the fractional abundance (Patterson & Fishbein, 1989). A sample was considered statistically significant if the total counts in the sample were greater than the probable error. All 30 samples were deemed statistically significant. The following formula was used to calculate standard error (S_Xi) for each species: ${\displaystyle S_{Xi}=1.96\sqrt{\frac{X_i\left(1-X_i\right)}{N_i}}}$ where X_i is fractional abundance of each species and N_i is the total species counts in a sample (Patterson & Fishbein, 1989). Species were included if their standard error was smaller than the fractional abundance. Based on these calculations, 32 species were deemed statistically significant and included in the calculations.

Species diversity was calculated using the Shannon Diversity Index (SDI) to explore changes in lake health over time (Shannon, 1948). The SDI is calculated using the formula: ${\displaystyle \mathrm{SDI}.=-\sum _{i=1}^S\left(\frac{X_i}{N_i}\right)\times ln\left(\frac{X_i}{N_i}\right)}$ where X_i is each taxon abundance in a sample, N_i is a sample’s total abundance, and S is the species richness of a sample. Samples were considered stable if the SDI was between 2.5 and 3.5, in transition if SDI was between 1.5 and 2.5, and stressed if SDI was below 1.5 (Magurran, 1988). Ratios between centropyxid and difflugiid species were calculated after the approach of Neville, McCarthy & MacKinnon (2010) (Table 1; Fig. S3). Stressed assemblages were defined as those with higher proportions of centropyxids to difflugiids, as centropyxid species have been found to be much more resistant to contaminated conditions as opposed to difflugiids, which have a low tolerance to contaminated conditions (Neville, McCarthy & MacKinnon, 2010; Nasser et al., 2016).

Geochemical data were screened prior to statistical analyses using the approach of Nasser et al. (2016) and Reimann et al. (2008). Variables that had missing values, or that were above or below the detection limit in more than 25% of their values were removed. All geochemical data was first converted to ppm and then normalized against aluminum, as it is a generally stratigraphically immobile element (Forghani et al., 2009). The LOI data were partitioned and reported as percentage water, organics, carbonates and minerogenics. Data reduction of the geochemical data was carried out using PRIMER (version 6) software by calculating Euclidean distance matrices for each element and running a 2Stage matrix (after Clarke & Gorley, 2005; Xu & Xu, 2016). Elements that were highly correlated and had no clear impact on arcellinidan assemblages were removed (after Fishbein & Patterson, 1993). Remaining elements were analyzed by using a RELATE function to run Spearman correlations, used for the typically skewed ecological datasets (i.e., datasets exhibiting a non-normal distribution), between elements and a Bray-Curtis similarity matrix of square root-transformed species counts (after Clarke & Gorley, 2005). Five variables (minerogenics, organics, arsenic, iron and mercury) that were found to be significantly correlated with species assemblages were used in the next stage of analyses (Table S2). Arsenic, mercury and iron were also compared on the basis of their ratio to sulfur and manganese to determine whether there was evidence of mobilization of elements resulting from redox reactions (Fig. 4). Further Spearman correlations were run between arsenic and elements it is often associated with; iron, manganese, and sulfur (Couture et al., 2013; Takamatsu, Kawashima & Koyama, 1985; Belzile & Tessier, 1990).

The make-up of the species assemblages was determined using CONISS cluster analysis—broken stick model in C2 (Grimm, 1987) (Fig. 5). Relationships between assemblages, core depth and arsenic were explored by creating a stratigraphic profile in C2 (Juggins, 2012) (Fig. 6). Non-metric multidimensional scaling (NMDS) function in RStudio (R Core Team, 2014) was employed to confirm the results of CONISS and assess the similarity between identified assemblages in multidimensional space (Kruskal, 1964). The one-way ANOSIM function in PRIMER was used to look for significant differences between assemblages (Clarke & Gorley, 2005). The five selected environmental variables were visually represented using line graphs created in Excel (Fig. 4).

Chronology and sedimentology

High-resolution chronological control indicated that all AMS ¹⁴C dates in the upper 15 cm were of material of dateable age but were highly jumbled and most likely of an allochthonous origin (e.g., old carbon bearing sediments removed from elsewhere and dumped into the lake; old carbon bearing sediments washed into the lake from the lake catchment). Dates are presented in years before present (¹⁴C age BP) in Table 2 (sample at 3 cm missing). The radiocarbon date at 19.5 cm–20.5 cm was accurate as it correlated well with the known stratigraphy from lower in the core. A ²¹⁰Pb profile from a nearby core in the same lake basin indicate that the upper 17 cm of the core dated from ∼1962 to the core collection date in 2012, and below an unconformity the underlying interval from 18 cm to 30 cm was determined to have been deposited during the early Holocene >7,000 years BP (Figs. 3 and 4; Table 3).

The 1FR core was characterized by darker sediments (5Y 2.5/1) in the top 44 cm of the core, with an ash layer occurring between 14 and 15 cm (Figs. 3 and 4). Further refining of the age model was possible due to the presence of three prominent stratigraphic markers comprised of wood ash with known depositional dates. These ash layers were derived from the burning of discarded Christmas trees in bonfires on the lake as part of New Year’s Eve celebrations in 1968, 1969 and 1971 (Figs. 3 and 4; Table 3) (V Sterenberg, former Yellowknife mayor, D Lovell, pers. comm., 2015). An unconformity at 17 cm marked the boundary between early Holocene sedimentation (7,230 ± 128) and sediments deposited after 1962 (Table 3; Figs. 3 and 4).

Loss on ignition

The loss on ignition (LOI) data (Table 1) showed an increase in both water and minerogenic content between the 6 and 17 cm core intervals. Organic content was correspondingly lower through this interval. Carbonate concentration remained stable throughout the entire 30 cm core interval studied. The minerogenic content were found to be the most significantly correlated environmental variable with arcellinida assemblages, with a Spearman’s Rho value of 0.317 and a p value of 0.002 (n = 30 samples) (Table S2; Fig. 4). Organics were calculated to be the second highest ranked environmental variable correlated with arcellinidan assemblages with a Spearman’s Rho value of 0.228 and a p value of 0.004 (Table S2; Fig. 4).

Arcellinida

(a) CONISS, ANOSIM and NMDS

The CONISS—cluster plot (Fig. 5) displayed three distinct assemblages at the 30% similarity mark. Assemblage 1 spanned a depth of 17–30 cm, assemblage two spanned a depth of 7–16 cm and assemblage three spanned a depth of 1–6 cm. Three distinct assemblages were identified by the NMDS bi-plot and matched CONISS results with the exception of sample FL5, with characteristics of both Assemblages 2 and 3, found at the stratigraphic boundary of these two assemblages. The results of the one-way ANOSIM test and comparison to the CONISS-defined species assemblages indicated that the assemblages were significantly different with global R statistic of 0.832 and a p value of 0.001 (n = 30 samples, n = 32 species). These results were further corroborated by the NMDS bi-plot, which revealed three groups of samples clustering closely to form the three identified arcellinidan assemblages (Fig. 7).

Geochemical variables

Arcellinida abundances were significantly correlated with arsenic (Table S2) with the third-ranked Rho value of 0.196 and p value of 0.034. Arsenic levels showed an increase over time, peaking at 15 and 8 cm, but showed a drop at 12 and 5 cm (Figs. 4 and 6). Arcellinida also showed a significant correlation with iron with an R-value of 0.192 and p value of 0.019 (Table S2). Iron levels peak at 24 cm but begin to show a marked decrease up-core of 20 cm (Fig. 4). Mercury was also found to be significantly correlated with arcellinidan abundances with an R-value of 0.168 and a p value of 0.046 (Table S2). Mercury values peak at 16 cm, then drop at 12 cm, but peak once again at 5 cm (Fig. 4). A comparison of mercury and arsenic distribution through the core to that of manganese, iron and sulfur indicated the second peaks of mercury and arsenic were possibly due to the mobilization of the elements (Fig. 4).

The correlation between arsenic and iron was significant (R = 0.336, p = 0.001, n = 140), as well as the correlation between arsenic and manganese (R = 0.093, p = 0.002, n = 140). The correlation between arsenic and sulfur was not found to be significant at R =  − 0.059, p = 0.954, n = 140.

Aerial photos

Aerial photography capturing eight decades of urbanization around Frame Lake indicated that while the area around the lake was relatively undisturbed in 1937, it had major development by 1964 (Fig. 2).

Assemblages

The results of the ANOSIM test, Spearman correlations, NMDS and CONISS cluster analyses demonstrate that Arcellinida sample assemblages could be statistically separated into three distinct groups, which correlated closely with their response to environmental conditions (Fig. 5; Table S2). Assemblage 1 (Baseline Limnological Conditions Assemblage [BLCA]) characterized samples deposited below the unconformity during the early Holocene. Assemblage 2 (Arsenic Contamination Assemblage [ACA]) was associated with mine waste influence and Assemblage 3 (Eutrophication Assemblage [EA]) developed during the more recent decades when deposition of mine contaminants ceased but level of nutrient input to the lake increased.

Baseline Limnological Conditions Assemblage—BLCA (n = 14)

The BLCA ranges from samples 17 to 30 and corresponds to ancient sediments deposited in the lake basin in the early Holocene, >7,000 years BP (Figs. 3 and 4). The SDI indicated that samples in this assemblage were all in transition, with the exception of samples 17, 18 and 28 that were in the stable range (Table 1) (Shannon, 1948; Patterson et al., 2013). Samples 17 and 18 mark the location of the unconformity between early Holocene and 20th century horizons in the lake, and sample 28 may mark an ecological transitional zone representing early hydroecological conditions as the lake transitioned from a basin within the very large Proto-Great Slave Lake into a small isolated lake. The difflugiid/centropyxid ratio revealed a strong dominance by difflugiids with a mean of 0.78 (Table 1; Fig. S3) (Neville, McCarthy & MacKinnon, 2010). Samples comprising the BLCA clustered together at both the 30% and 40% similarity marks, as seen in the CONISS plot and confirmed by the NMDS bi-plot (Figs. 5 and 7).

The BLCA is characterized by stable levels of minerogenics and organics, very similar to those found in the EA (Fig. 4, Table 1). Arsenic concentrations corresponding to the BLCA remain relatively stable, ranging from 289.8 to 121.4 ppm (Figs. 4 and 6). Mercury concentrations in the BLCA show a steady trend of short peaks and drops, but remain relatively constant throughout the interval (Fig. 4). Iron concentrations show a steady increase in the BLCA, peaking at 23 cm, but sharply dropping at 26 cm (Fig. 4). The faunal diversity of the BLCA is dominated by Difflugia protaeformis (Ehrenberg, 1830) strain “amphoralis” (22%), D. oblonga “oblonga” (19%) and D. oblonga “spinosa” (8%) (Table 1).

Arsenic Contamination Assemblage—ACA (n = 10)

The ACA spans samples 7–16 and corresponds to deposition between ∼AD 1962 and 1994 (Figs. 3 and 4; Table 3). The SDI indicates that these samples are in transition, with the exception of sample 16 that was in the stable range (2.6), marking an ecological transitional zone between the ACA and BLCA (Table 1) (Shannon, 1948; Patterson et al., 2013). The difflugiid/centropyxid ratio showed a significant drop in difflugiid abundances, with a mean value of 0.33 (Table 1; Fig. S3) (Neville, McCarthy & MacKinnon, 2010). The ACA clustered together in the CONISS plot (Fig. 5), indicating samples 7–16 are similar at both the 30% and 40% levels. Samples hosting this assemblage were also shown to exhibit tight grouping when examining the NMDS bi-plot. Interestingly, while sample FL7 appear to be part of this assemblage based on the CONISS results, the NMDS bi-plot revealed a notable break between the sample and the assemblage (Fig. 7). This could be attributed to the fact that the assemblage composition of FL7 is notably different than that of ACA.

The ACA is characterized by rising levels of minerogenics that peak at 11 cm and correspondingly low levels of organics that hit their lowest concentrations at 10 cm (Fig. 4, Table 1). Arsenic levels in the ACA show a double peak at 8.5 and 15 cm as well as a drop at 12 cm (Figs. 4 and 6). Mercury levels in the ACA show a double peak at each end of the assemblage interval and also have a marked drop at 12 cm (Fig. 4). Iron levels remain relatively constant throughout the ACA interval (Fig. 4). The faunal diversity of the ACA is characterized primarily by Centropyxis constricta (Ehrenberg, 1843) strain “aerophila” (23%), Centropyxis constricta (Ehrenberg, 1843) strain “constricta” Centropyxis aculeata (Ehrenberg, 1832), strain “aculeata” and Centropyxis aculeata (Ehrenberg, 1832), and strain “discoides” (Table 1).

Eutrophication Assemblage—EA (n = 6)

The EA is comprised of samples that range from 1 to 6 cm and corresponds to deposition in the lake from 1994 to 2012 (Figs. 3 and 4; Table 3). The SDI indicated that EA samples are in transition, with the exception of sample 5 that was in the stable range (Table 1) (Shannon, 1948; Patterson et al., 2013). The difflugiid/centropyxid ratio showed a dominance of difflugiid species with all values over 0.55 (Table 1; Fig. S3) (Neville, McCarthy & MacKinnon, 2010). The CONISS plot demonstrated that samples in the EA showed strong similarity as they clustered together at the 30% similarity range (Fig. 5). This similarity is also salient in the results of NMDS with the exception of sample FL5 plotting close to the ACA (Fig. 7). This is not surprising as the faunal composition of the sample is generally similar to that of the EA while also hosting a notable number of centropyxid species that are more characteristic of the ACA.

The EA corresponds to stable levels of organics, minerogenics and iron (Fig. 4, Table 1). Arsenic levels in the EA show a continuous decrease while mercury levels show a peak closer to the ACA and then drop further up core (Figs. 4 and 6). The faunal diversity of the EA is comprised mostly of C. tricuspis, C. constricta “aerophila”, D. oblonga “spinosa”, and D. oblonga “oblonga” (Table 1).