Historical phosphorus dynamics in Lake of the Woods (USA–Canada) — does legacy phosphorus still affect the southern basin?

ABSTRACT

Edlund MB, Schottler SP, Reavie ED, Engstrom DR, Baratono NG, Leavitt PR, Heathcote AJ, Wilson B, Paterson AM. 2017. Historical phosphorus dynamics in Lake of the Woods (USA–Canada) — does legacy phosphorus still affect the southern basin? Lake Reserv Manage. 33:386–402.

A historical phosphorus (P) budget was constructed for southern Lake of the Woods. Sediment cores from 7 bays were radioisotopically dated and analyzed for loss-on-ignition, P, Si, diatoms, and pigments. Geochemical records for cores were combined using focusing factors for whole-basin estimates of sediment, total P, and P fraction accumulation. Although historical monitoring shows that external P loads decreased since the 1950s, sediment P has continued to increase since the mid-20^th century. Much sediment P is labile and may be mobile within the sediments and/or available for internal loading and resuspension. Two mass-balance models were used to explore historical P loading scenarios and in-lake dynamics, a static one-box model and a dynamic multi-box model. The one-box model predicts presettlement external loads were slightly less than modern loads. The dynamic model shows that water-column P was higher in the 1950s–1970s than today, that the lake is sensitive to external loads because P losses from burial and outflow are high, and that the lake is moving to a new steady state with respect to water-column P and size of the active sediment P pool. The active sediment pool built up in the mid-20^th century has been depleted through outflow and burial, such that its legacy effects are now minimal. Comparison of historical nutrient dynamics and sediment records of algal production showed a counterintuitive increase in production after external P loads decreased, suggesting other drivers may now regulate modern limnoecology, including seasonality of P loading, shifting nutrient limitation, and climate warming.

KEYWORDS:

• Cyanobacteria
• internal phosphorus loading
• large lakes
• paleolimnology
• shallow lakes

Control and mitigation of excess nutrients, particularly phosphorus (P), continue to dominate lake management efforts (Schindler 2012, Schindler et al. 2016). In the USA, over 40% of lakes are impaired for phosphorus (USEPA 2016) and nutrient triggered cyanobacterial blooms are a global problem (Paerl et al. 2011). Measurements and models for determining basin P loading and sediment P burial, resuspension, and aerobic and diffusive loading are critical for addressing nutrient management and recovery from eutrophication. Many methods and models have been developed to estimate whole-basin and sediment fluxes (James and Barko 1993, Brenner et al. 2006), P retention capacity of sediments (Kopáček et al. 2007, Wilson et al. 2010), and long- and short-term P dynamics (Xie and Xie 2002, Wang et al. 2003, Norton et al. 2011). Importantly, these modeling exercises have been directed at nutrient-impaired waters throughout the world, although lake-specific models are often required (Havens et al. 2001). Resulting management efforts primarily target point and non-point P loadings; however, impaired lake conditions are often exacerbated by internal P loading through chemical release (especially under anoxic hypolimnetic conditions; Boström et al. 1982) and sediment resuspension. Internal loading may continue to determine lake condition even after significant reduction in external loads (Jeppesen et al. 2005, McCrackin et al. 2016).

Lake of the Woods (LOW) is a large, multibasin lake located along the borders of Minnesota (USA), Ontario and Manitoba (Canada). The lake extends about 100 km on both longitudinal and latitudinal axes, with the largest surface area in Big Traverse Bay, which connects to several secondary basins including Buffalo and Muskeg bays to the south and west, and Sabaskong Bay to the east. Water flows northward through Little Traverse Bay before passing through Big Narrows to join outflow from several deeper Canadian basins, bays, and outflows before discharging into the Winnipeg River at Kenora, Ontario. Overall, mean residence time is 1.71 yr (2000–2010; Zhang et al. 2013). Its major inflow is the Rainy River, which enters the southeast end of the lake near Baudette, Minnesota (Anderson et al. 2017).

With the publication of the Lake of the Woods State of the Basin Report (DeSellas et al. 2009; updated 2^nd edition: Clark et al. 2014) and the Minnesota Pollution Control Agency's placement of the lake in 2008 on the state's list of waters impaired for nutrients and eutrophication indicators, the future of the lake became a high-profile concern for Canada, Minnesota, First Nations and tribal governments, as well as the lake's stakeholders. The Basin Report highlighted nutrients and their biological impacts—primarily cyanobacterial blooms and a perceived increase in the frequency and extent of these nuisance blooms—as a primary resource concern for the lake.

Lake of the Woods has elevated concentrations of P in comparison to other lakes on the Precambrian Shield, a strong N–S gradient of water quality (Pla et al. 2005), and extensive cyanobacterial blooms (Binding et al. 2011). Although these characteristics have some historical precedence (Anderson et al. 2017), recent trends in lake ecology have been at odds with known effects of resource management. For example, monitored P loads from the Rainy River, the primary external source of P, have decreased over the last 40 yr, mainly due to improved management of point sources (Hargan et al. 2011). Following nutrient abatement programs, Rainy River water quality between the 1990s and 2000s shows little change in nutrient content (Hargan et al. 2011), which is further reflected in minimal change in in-lake nutrient concentrations based on the limited monitoring data available. Furthermore, paleoecological evidence from Canadian waters of northern LOW demonstrates little change in diatom-inferred P values (Rühland et al. 2008, 2010, Hyatt et al. 2011, Paterson et al. 2017), whereas fossil reconstructions from a small bay in the south basin show increasingly eutrophic conditions (Reavie and Baratono 2007). Cyanobacterial blooms are perceived to be more frequent and of greater spatial coverage than in previous decades, particularly in the southern basin, although evidence from monitoring, including satellite imagery, is equivocal (Chen et al. 2007, 2009, Binding et al. 2011).

Weak relationships between documented declines in nutrient influx and observed water quality may reflect either a strong legacy effect of sedimentary nutrients or establishment of alternative mechanisms regulating limnological conditions, such as climate-induced reduction in water-column mixing and reduced thermal structure (Paerl and Huisman 2008). In response to these challenges, management initiatives include an increase in the spatial and temporal resolution of monitoring, evaluation of satellite imagery, tests for cyanobacterial toxins, and development of a comprehensive P mass budget for the lake (Clark et al. 2014). To complement these initiatives, managers also need a detailed historical evaluation of nutrient dynamics of LOW to quantify the magnitude and timing of disconnect between changes in nutrient loading and lake response. In particular, sediments record changes in sedimentary P accumulation, as well as the chemical nature of P fractions, and often reveal how these factors vary in response to external loading, land use, climate, and other factors (Anderson et al. 1993, Ginn et al. 2012).

Historical and paleoecological techniques for estimating past P influx and temporal dynamics have proven useful in developing cooperative management strategies in other nutrient-enriched transboundary waters such as Lake Pepin and Lake St. Croix, smaller basins within the Upper Mississippi River (Edlund et al. 2009, Engstrom et al. 2009, Triplett et al. 2009). In those lakes, historical phosphorus mass balances, which estimated inputs based on the sum of whole-basin burial estimates and diatom-based estimates of P loss in outflows, indicated that P loading to each lake had increased rapidly after World War II in response to growing populations and increased point and nonpoint source loadings. Concomitantly, diatoms showed dramatic ecological changes in the last 200 yr, while diatom-inferred P concentrations increased after Euro-American settlement and the mid-20^th century. In contrast, recent analyses of sedimentary P, diatoms, and fossil pigments from phytoplankton in larger prairie basins (e.g., Lake Winnipeg, Lake Manitoba) suggest that lake production can be disconnected from estimates of P influx, particularly in poorly stratified waters (Bunting et al. 2016). Given the size and depth of southern LOW, it may be difficult to predict how production may have responded to nutrient management.

This project uses a combination of sedimentary P analysis, multi-proxy fossil analysis of phytoplankton (diatoms, pigments), and dynamic nutrient modeling to reconstruct historical changes in nutrient fluxes and conditions in southern LOW. In conjunction with a coupled paleolimnology effort (Reavie et al. 2017), we address these research questions:

1.	Does the sediment P record accurately reflect the lake's P loading history?
2.	How have P loadings to LOW changed over the last 150 yr?
3.	Can in-lake P dynamics be modeled to understand historical, legacy, and future nutrient dynamics?
4.	Do trends in core biogeochemistry and biological indicators reflect historical nutrient dynamics?

Materials and methods

Coring

Sediment cores were recovered from deep, flat depositional zones in 7 basins in LOW (Table 1; see also Reavie et al. 2017). Most cores were recovered with a piston corer consisting of a 6.5 cm diameter polycarbonate tube outfitted with a piston and operated with rigid drive rods from the ice surface (Wright 1991). A follow-up core was recovered the following summer from Buffalo Bay using a gravity corer (Renberg and Hansson 2008). Piston cores ranged in length from 90 to 98 cm, and the gravity core from Buffalo Bay was 9.5 cm long. All cores were stabilized with Zorbitrol (Tomkins et al. 2008) or sectioned immediately in the field in 0.5 cm increments to 10 cm depth using a vertical extrusion system. For piston cores, unextruded core material was sealed in its polycarbonate tube and transported horizontally back to the laboratory for further sectioning in 1 cm increments from 10 cm to 35 cm (to 60 cm for Sabaskong and Big Narrows cores).

Table 1. Lake of the Woods core names, dates, coring locations, depth at core site, and core recovery. Focusing factors are estimated by the flux of unsupported ²¹⁰Pb to the core site relative to known atmospheric depositional rates in the region (∼0.45 pCi/cm²/yr).

	Date	Lat (N)	Long (W)			Depth	Recovery	Focus
Core name	yyyymmdd	°N	°W	State/Province	Type	(m)	(m)	factor
LoW_BigNarrows	20120228	49.39472°	94.79395°	Ontario	Piston	8.53	0.98	1.68
LoW_LittleTrav	20120228	49.24643°	94.67145°	Ontario	Piston	9.18	0.98	1.37
LoW_Sabaskong	20120229	49.10064°	94.42108°	Ontario	Piston	6.85	0.98	1.87
LoW_BigTrav3	20120229	49.01931°	94.75391°	Minnesota	Piston	10.2	0.9	1.21
LoW_BigTrav4	20120301	49.08941°	94.99497°	Minnesota	Piston	10.13	0.96	1.44
LoW_Muskeg	20120301	48.97849°	95.17970°	Minnesota	Piston	8.08	0.95	1.59
LoW_BB2 H	20120818	49.10960°	95.22796°	Manitoba	HTH	5.52	0.095	0.41

Isotopic dating, biogeochemistry, and whole-basin deposition

Sediment cores were analyzed for ²¹⁰Pb activity to determine age and sediment accumulation rates over the past 150–200 yr. Lead-210 activity was measured from its daughter product, ²¹⁰Po, which is considered to be in secular equilibrium with the parent isotope. Aliquots of freeze-dried sediment were spiked with a known quantity of ²⁰⁹Po as an internal yield tracer and the isotopes distilled at 550 C after treatment with concentrated HCl. Polonium isotopes were then directly plated onto silver planchets from a 0.5 N HCl solution. Activity was measured for 1–3 × 10⁵ s using an Ortec alpha spectrometry system. Supported ²¹⁰Pb was estimated by mean activity in the lowest core samples and subtracted from upcore activity to calculate unsupported ²¹⁰Pb. Core dates and sedimentation rates were calculated using the constant rate of supply model (Appleby and Oldfield 1978, Appleby 2001). Dating and sedimentation errors represented first-order propagation of counting uncertainty (Binford 1990). For cores with problematic activity profiles, gamma spectrometry was used to measure supported ²¹⁰Pb (as ²¹⁴Pb) and identify the 1963 dating marker associated with the peak in ¹³⁷Cs activity. The short-lived isotope ⁷Be (half life 53.2 d) was also measured in the uppermost intervals of select cores using gamma spectrometry to determine the extent of sediment mixing from bioturbation and resuspension.

To understand whole-basin depositional rates for various constituents including dry bulk sediment and P fractions, a “focusing factor” was calculated for each core using the method of Engstrom and Rose (2013) and Hobbs et al. (2013) to normalize for downcore fluxes among basins. Focus factors estimate the degree to which each core site integrates sediment within a basin by comparing atmospheric flux to unsupported ²¹⁰Pb inventory at the core site. Atmospheric flux of ²¹⁰Pb in northern Minnesota is estimated at 0.45 pCi/cm²/yr (Lamborg et al. 2013). Sedimentation rates for individual basins were corrected for sediment focusing, the data for all cores pooled, and averaged among time intervals represented by approximately equal numbers of observations (5 yr window back to 1990, decadal intervals to 1940, 20 yr intervals to 1900, and pre-1900 samples grouped) to estimate whole-lake sedimentation rates.

Bulk-density (dry mass per volume of fresh sediment), organic, carbonate, and mineral content, and biogenic silica (BSi) concentrations and accumulation rates were determined for all cores. Details of these geochemical procedures are provided by Reavie et al. (2017). Sediment P was analyzed following the sequential extraction procedures in Engstrom (2005) and Engstrom and Wright (1984). Extracts were measured colorimetrically on a Lachat QuikChem 8000 flow injection autoanalyzer. Sediment P concentrations were also converted to flux using bulk sedimentation rates in each core. In addition to total P, sediment fractions include the refractory forms HCl-P and Organic-P and labile forms NaOH-P and exchangeable P (Ex-P).

Biological constituents measured in all cores included diatom and chrysophyte microfossils and fossil algal pigments; analytical procedures and results are presented by Reavie et al. (2017). To estimate historical water-column total P, or diatom-inferred total P (DI-TP), a diatom calibration set constructed by Hyatt et al. (2011) was applied to relative abundance data of downcore diatom assemblages using weighted averaging regression with inverse deshrinking. Calibration model performance and reconstruction statistics are presented in Reavie et al. (2017).

Modeling historical phosphorus dynamics

Two modeling approaches were developed and applied to downcore data to understand historical nutrient dynamics, historical P loads, and current nutrient trajectories. Model 1 is a simple one-box whole-lake mass balance, whereas Model 2 is a 3-box dynamic model run from 1850 to present. Each model is presented below with its conceptual basis, assumptions, input data, and a discussion of its results, trends, potential shortcomings, and key findings. Model 2 was assembled and run using the software Stella 9.0 (isee systems, Lebanon, NH, USA; www.iseesystems.com).

Supporting data for modeling of historical P budgets came from several sources. For Model 1, historical water-column P was estimated using diatom-inferred total P (DI-TP) reconstructions from all cores (Table 2; Reavie et al. 2017). Lake area was calculated from polygons digitized from aerial photography using the software QGIS (QGIS Development Team 2015) and lake volume by basin was taken from Zhang et al. (2013). Outflow rates to the Winnipeg River at Kenora were available from 1927 to 2008 and provided by the Lake of the Woods Water Control Board (www.lwcb.ca). Outflow at the Big Narrows was scaled based on supplemental data provided in Zhang et al. (2013) by comparing daily step outflow from 2000 to 2010 at the Big Narrows to Kenora. Phosphorus loadings from the Rainy River were assembled from available records from 1954 to present, including compilations by Beak Consultants Ltd. (1990) and Hargan et al. (2011), and recent monitoring coordinated by the Minnesota Pollution Control Agency (Table 3). Data were summarized using decadal average flows and arithmetic means of measured TP. Other sources of P loads to the lake including atmospheric deposition, minor tributaries, and shoreline erosion were taken from Hargan et al. (2011) and HEI (2013).

Table 2. Model 1 output where I = B + O, P Inputs (I), P Burial (B), and P Outflow (O) are in tonnes P/yr. P Outflow is estimated from diatom-inferred TP (Reavie et al. 2017) multiplied by outflow volume (see Table 3).

Time interval	P input	P outflow	P burial
(years)	(t P/yr)	(t P/yr)	(t P/yr)
2005–2011	1326	361.0	965
2000–2004	1080	270.4	809
1995–1999	1000	278.5	721
1990–1994	908	243.7	664
1980–1989	806	219.7	586
1970–1979	753	286.1	467
1960–1969	701	279.6	422
1950–1959	935	221.2	714
1940–1949	859	215.1	644
1920–1939	712	147.2	564
1900–1919	676	153.0	523
pre-1900	559	128.2	431

Table 3. Model 2 input data, parameters, and data sources.

Parameter	Value	Units	Sourcea
Surface area of lake	2.83 × 10^9	m^2	GIS
Volume of lake	18.48 × 10^9	m^3	1
Outflow at Big Narrows
1850–1930	10.7 × 10^9	m^3/yr	1, 4
1940	13.7 × 10^9	m^3/yr	1, 4
1950	12.9 × 10^9	m^3/yr	1, 4
1960	15.8 × 10^9	m^3/yr	1, 4
1970	15.6 × 10^9	m^3/yr	1, 4
1980	12.8 × 10^9	m^3/yr	1, 4
1990	14.4 × 10^9	m^3/yr	1, 4
2000–2050	14.9 × 10^9	m^3/yr	1, 4
Phosphorus load from Rainy River
1850–1900	300	t/yr	5
1900	400	t/yr	5
1910	500	t/yr	5
1920	600	t/yr	5
1930	800	t/yr	5
1940	1000	t/yr	5
1950	1176	t/yr	3
1960	1319	t/yr	2, 3
1970	830	t/yr	2, 3
1980	546	t/yr	2, 3
1990	519	t/yr	2, 3
2000	377	t/yr	2, 3
2000–2050	350	t/yr	2, 3
Phosphorus load from other sources
1850–2050	232	t/yr	2, 6
Whole-basin sediment accumulation rate (areal)
1850–1900	0.238	kg/m^2/yr	5
1900–1919	0.288	kg/m^2/yr	5
1920–1939	0.335	kg/m^2/yr	5
1940–1949	0.318	kg/m^2/yr	5
1950–1959	0.341	kg/m^2/yr	5
1960–1969	0.353	kg/m^2/yr	5
1970–1979	0.384	kg/m^2/yr	5
1980–1989	0.417	kg/m^2/yr	5
1990–1999	0.469	kg/m^2/yr	5
2000–2050	0.506	kg/m^2/yr	5

^a Supporting data: (1(1) $$\mathrm{I}=\mathrm{B}+\mathrm{O}$$(1) ) Zhang et al. 2013, (2(2) $$\mathrm{Cumulative}\mathrm{P}\mathrm{in}\mathrm{active}\mathrm{layer}=\mathrm{Ext}\mathrm{Load}\times \%\mathrm{to}\mathrm{Sed}-\mathrm{Burial}-\mathrm{InLoad}$$(2) ) Hargan et al. 2011, (3(3) $$\mathrm{Cumulative}\mathrm{buried}\mathrm{P}=\left(\mathrm{Cum}.\mathrm{P}\mathrm{in}\mathrm{active}\mathrm{layer}/\mathrm{MS}\right)\times \mathrm{Sed}\mathrm{Rate}$$(3) ) Beak Consultants Ltd. 1990, (4(4) $$\mathrm{Lake}\mathrm{P}=\mathrm{Ext}\mathrm{Load}\times \left(1-\%\mathrm{to}\mathrm{Sed}\right)+\mathrm{InLoad}-\mathrm{Out}$$(4) ) Matt DeWolfe (lwcb.ca), (5) this study, (6) HEI 2013.

Results

Sediment core records

Most cores from LOW showed monotonic exponential declines in ²¹⁰Pb inventories to depths with background (supported) activities (Fig. 1). Cores generally reached supported levels of ²¹⁰Pb around 25–35 cm depth, except for Buffalo Bay, where supported levels were reached at 7–8 cm. Supported ²¹⁰Pb activities ranged from 0.85 pCi/g (Muskeg Bay) to 1.28 pCi/g (Big Narrows). Sediments dated to 1900 correspond to the approximate period of European settlement and damming of the lake at Kenora (Clark et al. 2014) and were found between 17 cm (Little Traverse) and 34 cm (Sabaskong Bay) downcore, except for Buffalo Bay (∼7.5 cm). Buffalo Bay began to accumulate lacustrine sediments at ca. 1900, likely in response to damming at Kenora, which raised LOW water levels by ∼1 m (Clark et al. 2014). Sediment focusing factors varied among the core sites from 0.41 at Buffalo Bay to 1.87 in Sabaskong Bay (Table 1). The short-lived isotope ⁷Be (half-life 53.2 d) was measured in select cores and detected to depth of 1–4 cm; if ⁷Be can be detected in sediments dated by ²¹⁰Pb at 6–10 yr old, sediment mixing must be occurring in LOW at least to some degree (data not shown).

Figure 1. Downcore profiles for 7 Lake of the Woods cores for total ²¹⁰Pb activity, date–depth relationship, and sedimentation rate plotted against core depth (cm). Dashed line in ²¹⁰Pb inventory represents level of supported ²¹⁰Pb.

Most cores showed increasing sedimentation rates in more recent deposits with modern rates typically 2-fold greater than those before 1900 (Fig. 1). Some cores had slightly greater increases in sedimentation rates including the Big Traverse Bay and Little Traverse Bay cores, with recent sedimentation nearly 3 times pre-1900 rates. Little Traverse and Muskeg bays had secondary increases in sedimentation rates since the 1970s and 1980s, respectively. Modern sedimentation rates varied from 0.6 (Big Traverse 4) to 1.2 kg/m²/yr (Sabaskong Bay), whereas presettlement rates ranged from less than 0.1 (Big Traverse 3) to 0.6 kg/m²/yr in Muskeg Bay. Following correction for sediment focusing in each basin and pooling of all cores based on averaged time intervals, estimates of whole-lake sedimentation rates increased from a presettlement rate of 0.27 kg/m²/yr to a peak in the 1970s of 0.69 kg/m²/yr. Whole-basin sedimentation rates declined slightly in the 1980s but have risen to approximately 0.7 kg/m²/yr since the 1980s (Fig. 2a).

Figure 2. Whole-basin estimates of focus corrected sediment accumulation and diatom-inferred historical water column total P plotted against time period. Fig. 2a. Whole-basin estimates of focus corrected sediment accumulation (kg/m²/yr). Fig. 2b. Whole-basin estimates of water column diatom-inferred total P (DI-TP; µg/L).

Total P in LOW sediment ranged from 0.4 to over 1.0 mg P/g dry mass (Fig. 3). The organic-P and NaOH-P fractions were most abundant in Big Traverse 4, Little Traverse Bay, Sabaskong Bay, and Big Narrows. In contrast, HCl-P was a predominant P fraction in Big Traverse 3, Buffalo Bay, and deeper sediments of Little Traverse and Muskeg bays. In all cores the accumulation rates of sediment P and fractions increased 2- to 3-fold over the 20^th century, with the highest levels at the core surface. Based on historical estimates of P loading from the Rainy River, there have been significant declines in P loading since the mid-1970s to present day that are 2- to 3-fold less than loading estimates derived from 1950s–1970s. However, there is no clear indication of decreased accumulation of P in the sediments in response to decreased external loads, possibly because upward mobility of P within the sediments obscures the trend of P inputs to the sediments (James et al. 2015).

Figure 3. Geochemistry of 7 Lake of the Woods cores including concentration (mg P/g sediment) and flux (mg P/cm²/yr) of total sediment phosphorus and phosphorus fractions including HCl-P, NaOH-P, Organic-P, and Exchangeable-P, and water column diatom-inferred total phosphorus (DI-TP; µg/L) estimates from Reavie et al. (2017) plotted against core date.

Whole lake P accumulation rates were estimated from the time-averaged sum of P accumulation estimates from all sites, each independently corrected for sediment focusing (Fig. 4a). The P fractions were also treated separately as refractory (HCl-P, Org-P) or labile (potentially exchangeable) fractions (Ex-P, NaOH-P; Fig. 4a). Labile fractions are prevalent in all levels of LOW sediments with the amount increasing upcore, consistent with the expectation that these P fractions are potentially mobile within the sediment profile.

Figure 4. Whole-basin estimates of historical accumulation of phosphorus (P) and P fractions in Lake of the Woods sediments by time period. Fig. 4a. Accumulation of P differentiated into refractory components (HCl-P and Organic-P; green bars) and labile components (NaOH-P and Exchangeable-P; yellow bars); minimum burial estimates of refractory fractions were used in Model 2. Fig. 4b. Conceptual model of the Active and Buried inventory of P present in 2011 (see text for details).

Because burial of P is often the primary mechanism that removes P from a lake, we developed a conceptual model that considers the historically or permanently buried P and the active pool of P (Fig. 4b). We recognize that a significant proportion of the P in upper sediment layers represents an active pool of P that can be exchanged with the overlying waters or within the cores via mobility and bioturbation. In addition, the active pool is not restricted to the labile fractions because of resuspension (James 2017) and because labile P fractions are present in deep sediments (Fig. 4a). We also recognize from the ⁷Be inventory that sediments may be rapidly mixed in LOW down to 5 cm. Because of these factors (mixing, resuspension, within-core mobility) we do not know at the time of coring and at a given sediment depth what proportion of P is actually buried. Therefore, for modeling purposes our conceptual basis recognizes that there is a pool of P available for exchange (“Active”; Fig. 4b) and a pool of P that is truly buried and no longer available for exchange with the lake (“Buried”; Fig. 4b). Model 2 explores the behavior of these pools, particularly the net flux of P from the active pool via diffusion and resuspension to estimate water-column TP concentration, and uses the whole-basin inventory of sediment P (active plus buried) in sediments deposited from 1860 to 2011 as a modeling target (see Model 2 below).

Among the 7 cores analyzed for diatoms, most show continuous upcore increases in DI-TP (Fig. 3; see also Reavie et al. 2017). Analysis of all cores, except Buffalo Bay (no 19^th century sediments), suggested that background (pre–Euro-American settlement) TP concentrations in the water column to be approximately 10 µg P/L throughout the southern LOW. Cores from Muskeg, Big Narrows, and Big Traverse 4 showed increasing DI-TP upcore after 1900, whereas Big Traverse 3, Sabaskong, and Little Traverse had more marked increases in DI-TP after 1950. Overall, Buffalo Bay had the highest DI-TP values than all other cores from LOW with recent values exceeding 30 µg P/L. Values of DI-TP from the most recent sediments of other cores were typically between 20 and 30 µg P/L with several cores exceeding 30 µg P/L in the uppermost sections (Big Narrows, Muskeg, Buffalo Bay).

The DI-TP reconstructions of 6 cores were combined (Buffalo Bay omitted in pre-1900 as it did not preserve a pre-damming record) by time increment to estimate whole-lake historical water-column TP (Fig. 2b). Whole-lake DI-TP trends suggest TP concentration was about 10 µg P/L, which steadily increased to a peak of ∼18 µg P/L in the 1970s. The DI-TP estimates appear to be low compared to available monitoring data from the late 1960s, which indicate south basin TP concentrations of 30–100 µg P/L (Reavie et al. 2017). After the 1970s, DI-TP values remained between 15 and 17 µg P/L until the most recent period (2005–2011) when whole-lake DI-TP increased to over 24 µg P/L. Comparison of DI-TP with monitored TP values from within the cored basins suggest that average TP from 2005 to 2011 was 38 µg P/L and 31 µg/L in 1999 based on roughly monthly late spring–summer sampling during focused monitoring efforts by US and Canadian agencies. It is also apparent from the monitoring data that in the southern basins there were distinctly higher TP readings in the late summer months (>40 µg P/L) compared to spring (20–32 µg P/L; Lake of the Woods Water Sustainability Foundation 2011, Reavie et al. 2017). Whole-lake DI-TP (or for Model 2, calculated P concentration) was multiplied by discharge at Big Narrows, which was estimated from 1900–2011 based on scaling daily step outflows taken at both Kenora and Big Narrows from 2000–2010 (Fig. 5a; Zhang et al. 2013).

Modeling historical P dynamics

Two whole-basin modeling approaches were used to explore historical P loading scenarios to LOW and in-lake nutrient dynamics.

Model 1) Simple whole-lake mass balance

We first applied a commonly used one-box whole-lake mass flux model to estimate historical P loading in LOW (Rippey and Anderson 1996, Engstrom et al. 2009, Triplett et al. 2009, Engstrom and Rose 2013): (1) $$\mathrm{I}=\mathrm{B}+\mathrm{O}$$(1) where all external inputs (I) of P are either permanently buried in sediments (B) or removed from the lake via outflow (O). The sum of burial and outflow at any time is a first-order estimate of historical P loading to the lake. Modeled outflow (O) is estimated using the whole-lake historical diatom-inferred concentrations of TP (DI-TP; Fig. 2b) multiplied by the outflow at Big Narrows (Fig. 5). Whole-lake burial (B) of P was calculated from focus-corrected flux rates of total sediment P for each sub-basin as above (Fig. 4a). Burial of P is assumed to be permanent with only minor internal loading and no mobility within sediments (i.e., observed sediment flux reflects actual burial rate at each dated interval).

Model 1 P loading estimates for LOW are estimated to be approximately 559 t P/yr before settlement (Table 2). Modern whole-lake load estimates (based on monitoring) are only slightly higher and range from 582 t P/yr (2005–2014; RESPEC, unpublished) to 687 t P/yr (2005–2011; Hargan et al. 2011). After settlement, model results suggest P loadings increase continuously to modern rates of 1326 t P/yr (Table 2). Based on monitored loading estimates (see Hargan et al. 2011, Anderson et al. 2013, Zhang et al. 2013), this model clearly overestimates modern loadings to the lake. Importantly we also do not see any modeled decrease in loadings to the lake since the 1980s that would reflect well-documented decreases in P inputs from the Rainy River (Hargan et al. 2011). A large overestimate of modern P loads to the lake and no indication of decreased loading after 1980 (Fig. 4a) reflect shortcomings of this model and limit its applicability to sediment records deposited during steady state conditions during presettlement times. The assumption that LOW rapidly and permanently removes external P from the lake via burial is likely violated due to the within-core P mobility, high rates of resuspension, and slow sedimentation rates.

Model 2) Dynamic 3-box model with annual time step, 1860–2011

To better estimate temporal changes in TP influx and in-lake fluxes, a 3-box dynamic model was constructed and run from 1860 to present (Fig. 6). In this case, modeled pools (inventories) of P include buried sediment P (Cumulative buried P, Fig. 4b), an active sediment pool of P (Cumulative P in active layer, Fig. 4b) available for exchange with the water column or burial, and P in the water column (Lake P) from external and internal loading that are estimated using: (2) $$\mathrm{Cumulative}\mathrm{P}\mathrm{in}\mathrm{active}\mathrm{layer}=\mathrm{Ext}\mathrm{Load}\times \%\mathrm{to}\mathrm{Sed}-\mathrm{Burial}-\mathrm{InLoad}$$(2) (3) $$\mathrm{Cumulative}\mathrm{buried}\mathrm{P}=\left(\mathrm{Cum}.\mathrm{P}\mathrm{in}\mathrm{active}\mathrm{layer}/\mathrm{MS}\right)\times \mathrm{Sed}\mathrm{Rate}$$(3) (4) $$\mathrm{Lake}\mathrm{P}=\mathrm{Ext}\mathrm{Load}\times \left(1-\%\mathrm{to}\mathrm{Sed}\right)+\mathrm{InLoad}-\mathrm{Out}$$(4)

Input data for Model 2 are the external P loads (Ext Load) from the Rainy River, which were estimated annually from the 1950s to 2011 (Table 3), and other sources of P (other tributaries, shoreline erosion, atmospheric deposition), which were held constant from 1850 to 2011 at 232 t P/yr (Table 3). Initial external load conditions (1850–1900) were set at 300 t P/yr from the Rainy River plus 232 t P/yr from other sources (total external load 532 t P/yr), similar to Model 1 presettlement loading estimates (Table 3). From 1900 to 1950, P loads were increased incrementally to 1950s monitoring estimates (Table 3). The model also incorporated a 10 yr lag in burial; P that reached the sediments could not be permanently buried for 10 yr, but remained available for exchange with the water column as supported by the depth of mixing of ⁷Be and data from other large lakes (Nürnberg and LaZerte 2016).

Figure 5. Outflow and P loss at Big Narrows. Fig. 5a. Historical flows at Big Narrows for each time period, km³/yr. Fig. 5b. Estimates of historical loss of phosphorus through outflow at Big Narrows by time period (lower panel). P loss represents the whole-lake historical diatom-inferred total phosphorus multiplied by historical flows at Big Narrows for each time period.

Figure 6. Model 2 is a 3-box dynamic model run from 1850 to 2050. Three inventories of P are estimated including P in the lake (Lake P), Cumulative P in the Active Layer, and Cumulative P in the Buried Layer by adjusting the percent of external P load (EX) that goes to the sediment (% to Sed), the internal load rate (InLoad), and the mass of sediment (MS) that is in the Active Layer.

Model variables that were manipulated included the percent of external load that goes directly to sediment (% to Sed), which ranged from 0 to 50%, based on our knowledge that much of the P load from the Rainy River is in dissolved forms and readily available for in-lake production. The mass of sediment in the active layer (MS), or the mass of sediment in the top 0–10 cm depth increment, which ranged from 8.03 to 19.23 kg/m². MS represents the amount of sediment in the layer that can exchange P with the lake before becoming buried. The mass of sediment and P in this active layer determines the concentration of P at the time of permanent burial. The internal loading rate (InLoad) was also manipulated and represents a net annual flux calculated as the % of P in the active layer that enters the lake through resuspension and/or redox cycling and diffusion; range 0–2.5%.

Model variablesFigure 6., Figure 7.were manipulated through trial and error to best meet model target criteria (Table 3). First, the model was evaluated against known or modeled in-lake concentration of TP with targets set at 10 µg P/L presettlement based on whole-basin DI-TP (Reavie et al. 2017), 1960s TP monitored at approximately 70 µg P/L, and 2005–2011 TP values using whole-basin DI-TP of 25 µg P/L (Reavie et al. 2017). The second model target was the whole-basin inventory of P measured in sediments of southern LOW deposited in sediments from 1860 to 2011 (106,620 t P) and 1940 to 2011 (67,746 t P).

Target criteria were best satisfied when (1(1) $$\mathrm{I}=\mathrm{B}+\mathrm{O}$$(1) ) % to Sed was set at 75%, a reasonable number given that at least a quarter of TP entering LOW from the Rainy River is dissolved P; (2(2) $$\mathrm{Cumulative}\mathrm{P}\mathrm{in}\mathrm{active}\mathrm{layer}=\mathrm{Ext}\mathrm{Load}\times \%\mathrm{to}\mathrm{Sed}-\mathrm{Burial}-\mathrm{InLoad}$$(2) ) the InLoad was set at 2.5% of the Active Pool of P; and (3(3) $$\mathrm{Cumulative}\mathrm{buried}\mathrm{P}=\left(\mathrm{Cum}.\mathrm{P}\mathrm{in}\mathrm{active}\mathrm{layer}/\mathrm{MS}\right)\times \mathrm{Sed}\mathrm{Rate}$$(3) ) the active layer was defined as the top 0–5 cm of the core with a corresponding sediment mass (MS) of 8.03 kg/m². Model 2 results are presented from 1860 to 2011 (the model was run from 1850 to 1860 to reach initial steady state conditions, and extended to 2050 using current loading rates) and are best interpreted by examining model estimates of water-column TP and the size of the active pool of P (Fig. 7).

Figure 7. Output of Model 2, 1860–2050. Model 2B (blue line) is based on lower external inputs in presettlement times compared to Model 2A (red line; see text). Fig. 7a. Modeled water-column TP (µg/L) peaks in 1950–1960s with rapid water quality improvement after 1960s (blue line). Model 2B delays the rise of TP until 1900 (red line) but has no effect post-1960s. Stable TP levels are reached by 2015–2020 if modern external loads remain constant. Fig. 7b. Modeled P in the active layer of Lake of the Woods sediments. The active pool of P was greatest in the 1960s regardless of presettlement external load scenarios, and legacy P has been rapidly reduced since the 1960s. Model output suggests the active pool of sediment P is reaching a stable condition.

Model 2A overestimates initial water-column TP in LOW at just over 20 µg P/L, shows a rapid increase to peak levels of 77 µg P/L in the 1950s, and then depicts slowly decreasing TP to modern levels of 26 µg P/L. The active pool of P increases rapidly after 1900 to maximum levels in the 1960s before declining to modern levels by the 2010s. Two modifications were made to better understand model performance and future water quality trends. The model is highly sensitive to external loads. Hence, if external loads are reduced to 232 t P/yr from 1850 to 1900 (equivalent to current sources of P other than the Rainy River), Model 2B output more closely matches our presettlement DI-TP estimate of ∼10 µg P/L. Importantly, the increase in water-column TP is delayed until about 1900, which aligns with the timing of settlement, land use changes, and damming. If the model is run through 2050 holding external loads from the Rainy River at current estimates of 350 t P/yr, the lake reaches a steady state by 2020 with TP of 25 µg P/L and an active pool of 10,000 t P.

Overall, Model 2 shows water-column TP concentrations were 2–3 times greater in the 1950s–1970s than today, and that decreased external loading after the 1970s resulted in significant decreases of P concentration in the lake compared to the mid-20^th century. The lake is responsive to external loads because P burial and outflow are large net annual losses in LOW. Similarly, the active pool of sediment P was largest in the 1960s and that legacy pool of P has been rapidly depleted through burial or outflow to its current size of 10,000 t P. As such, the lake will approach a new steady state with regard to water-column TP and its active pool of P if current loading trends continue.

Discussion

Paleolimnological analysis of sediment cores is widely used in lake management to determine background or reference lake condition, periods and direction of lake change, an understanding of potential drivers of change, and current ecosystem trajectories (Smol 2009). In LOW, the paleolimnological approach was extended from a historical account of lake water quality and ecological consequences (Reavie et al. 2017) to a whole-lake interpretation of the stratigraphy of sediment P to more fully understand historical patterns of nutrient loading, quantify temporal variability in lake-sediment P dynamics, and evaluate current trends in lake conditions using traditional and dynamic modeling techniques. We organize our discussion of core records and modeling results based on our initial research questions followed by the limitations and management implications of this approach.

Does the sediment P record accurately reflect the lake's P loading history?

Historical observations suggest that TP influx to LOW has declined from maxima during the mid-20^th century. For example, estimates of TP influx compiled by Beak Consultants Ltd. (1990) and Hargan et al. (2011) rigorously account for monitored TP loads from the Rainy River as well as other tributary loads and sources during 1954–2011 (Table 3). These data suggest that Rainy River P influx was greatest during the 1950s (∼1700 t P/yr) but dropped by the 1970s, with a steady decline to modern loadings that range from 237 to 559 t P/yr (Table 3; Zhang et al. 2013). At the same time, P from smaller tributaries, atmospheric deposition, and shoreline erosion accounts for an additional 232 t P/yr (Hargan et al. 2011, HEI 2013).

Sediment P profiles in LOW do not directly record the dynamic nature of P influxes since ca. 1950. Instead geochemical analyses show the burden of P retained in the sediment is mobile. Its gradual upcore diffusion increases the amount of P observed in the upper sections of all cores and obscures the historical loading peak of the 1950s–1970s. This phenomenon is not uncommon in lake sediment cores from eutrophic lakes, especially those with relatively low sedimentation rates and with a higher propensity for recycling of sedimentary P into the water column (Carey and Rydin 2011, Ginn et al. 2012). In contrast, lakes with high sedimentation rates and rapid P burial can preserve known temporal patterns of historical P influx (Engstrom et al. 2009, Triplett et al. 2009), and cores will maintain that record based on repeat coring efforts separated by decades (Søndergaard et al. 2003, Blumentritt et al. 2013).

How have P loadings to LOW changed over the last 150 yr?

In LOW, a combination of paleolimnology, modeling, and monitoring was required to understand that P loadings were estimated to have increased rapidly following settlement to peak levels in the 1950s–1970s, after which loadings decreased rapidly following nutrient abatement regulations. Past changes in P influx in the absence of monitoring data have been estimated using a combination of whole-lake estimates of P burial and diatom-inferred estimates of water-column TP. For example, this approach has proven successful in developing nutrient and sediment reduction strategies in large transboundary lakes such as the Upper Mississippi River's Lake Pepin and Lake St. Croix (Edlund et al. 2009, Engstrom et al. 2009, Triplett et al. 2009). In these lakes, relatively high sedimentation rates provide rapid and efficient burial of P and a sediment record that reflects trends in P loading. However, because LOW sediments do not preserve a direct record of P loading, we cautiously applied a simple whole-lake mass balance model to estimate presettlement loadings to LOW. If we assume that the presettlement sediment record in LOW represents a long-term steady state, our Model 1 predicts presettlement P loading at 559 t P/yr. Because of upcore mobility of P in the sediments, Model 1 is limited in its application to presettlement (steady state) conditions. For other historical loading estimates, we must rely on monitoring data, which suggest peak loading from the Rainy River in the 1950s, slight declines through the 1970s, and a rapid decrease in loadings from the 1980s to present. Other modern sources of P are estimated at 232 t P/yr and include inputs from minor tributaries, atmospheric deposition, and shoreline erosion (Hargan et al. 2011, HEI 2013).

Can in-lake P dynamics be modeled to understand historical, legacy, and future nutrient dynamics?

Model 2 explored the historical behavior of P in LOW that led to the modern distribution of sediment P. This model was necessary because the abundance and distribution of P fractions in LOW sediment cores indicate there is a pool of readily exchangeable P, and that pool of P increases at the top of the core. This pattern was clearly identified in all cores in this study and by James (2017) from sites in Big Traverse and Muskeg bays. Because sediment P is potentially mobile, the amount of P at a particular depth (and therefore time) is transient. If a core is collected from LOW today, the downcore abundance of P is only a snapshot of current sediment P distribution, and that distribution reflects historical loading and in-lake processes that control P loading (internal and external), deposition, mobility, and burial. Likewise, a core taken in 1970 would have a different profile than today's core, and the interval dated from 1970 in today's core will not look like it did in 1970 in geochemical terms.

Whereas many modeling efforts strive to disentangle P dynamics at the sediment–water interface and within the oxic/anoxic sediment boundary (e.g., Wang et al. 2003), our model uniquely considered P dynamics at annual time steps on time frames greater than a century.

Model 2 results yield new insights on historical nutrient dynamics in LOW and provide perspective on current and future water quality trends in the lake. First, water-column P was significantly higher in the past, particularly in the 1950s–1970s, than it is today. Second, the lake is very responsive to changes in external loads. Model results show the lake quickly became more eutrophic as nutrient loading ramped up following settlement, but also show that water-column P levels quickly fell as external loads were reduced after the 1970s. No long-term trend in outflow volume and P loss at Kenora was noted that might account for this drop in water-column P (Table 3). Third, the responsiveness of the lake is a consequence of rapid and large burial and outflow fluxes that remove P from the lake. Last, with rapid reduction of external loads after the 1970s and current external loads remaining relatively constant for the last decade, LOW has both rapidly depleted any legacy pool of sediment P and has or will soon reach a new steady state with respect to water-column P and the size of its active pool of sediment P.

Do trends in core biogeochemistry and biological indicators reflect historical nutrient dynamics?

Biological remains preserved in the sediments of LOW record how ecological conditions changed in the lake over the last 150 yr in response to changing nutrient dynamics; however, the indicators of historical algal productivity in LOW sediments offer somewhat conflicting scenarios that need to be reconciled with our model reconstructions of historical P loading and dynamics. Community changes in the diatoms are presented in detail elsewhere (Reavie et al. 2017) and in conjunction with biogenic silica and fossil algal pigments provide a record of historical diatom productivity. Historical changes in cyanobacteria communities and productivity are similarly recorded by their fossilized pigments.

Pigment profiles, particularly those of general algal indicators (e.g., lutein-zeoxanthin) and diatom specific pigments (e.g., diatoxanthin) suggest 2 periods of high productivity in the recent history of LOW. The first period occurred from the 1950s through 1970s, during the peak of nutrient influx to LOW, and was followed by a decline in productivity in the 1980s followed by a second period of increased diatom productivity since the 1990s. There are significant changes in diatom communities in the most recent decades, particularly a greater abundance of species with higher TP optima including Cyclostephanos dubius, several small Stephanodiscus species, and Aulacoseira granulata (Reavie et al. 2017). This most recent diatom community represents a species assemblage not previously seen in the lake. Despite evidence from pigment proxies that suggest greater diatom productivity in the 1950s–1970s, there is no indication that the most recent high-P indicator taxa were common in the 1950s–1970s. As such the DI-TP does not effectively predict elevated P levels that were measured in the 1950s–1970s in LOW (Reavie et al. 2017: Supplement C). Similarly, biogenic silica records, whether treated as a concentration or flux, do not show increased diatom productivity during the 1950s to 1970s, even though external P loading to the lake was higher and diatom pigment indicators suggest higher productivity at that time (Reavie et al. 2017). Biogenic silica is normally treated as a proxy for historical diatom productivity, but in LOW produces a confounded record that is difficult to reconcile with sediment pigments and historical P loading.

Fossil pigments also indicate 2 periods of elevated cyanobacterial production in LOW. The first period is from the 1950s to 1970s and is characterized by high concentrations of cyanobacterial (e.g., echinone and canthaxanthin) and general algal indicators (e.g., lutein-zeoxanthin; Reavie et al. 2017). The same pigment groups show a second increase since the 1990s in most cores. However, there is also an increase since the 1990s of an additional pigment, myxoxanthophyll, an indicator of filamentous and colonial cyanobacteria including several of the potentially toxic forms (e.g., Microcystis), further suggesting that the biological communities present in the most recent decades are unique in the recent history of LOW.

Recent biological changes in LOW seem paradoxical in relation to the simple reduction of external P loads and depletion of the active pool of P as indicated by P monitoring and our modeling exercise. This incongruity suggests other factors must be driving changes in the algal communities. One potential driver is a shift in nutrient limitation. The few historical monitoring data on open-water nutrient stoichiometry suggest that the lake was P-limited in the 1960s and that reduction of point-source inputs has reduced N in a disproportionate ratio (relative to the Redfield ratio) to P leading to N-limitation (Pla et al. 2005, Reavie et al. 2017), an environmental factor linked to enhanced cyanobacterial production (Ferber et al. 2004, Orihel et al. 2012). Second, nutrient abatement efforts targeted point source loads (principally the pulp/paper industry and wastewater treatment plants), which has changed the seasonality of external loading to the lake from the Rainy River from more constant loading to maximum loading occurring April–June (J. Anderson, pers. comm.), likely affecting algal seasonality in the lake. Third, climate warming may have exacerbated gains in water quality made through nutrient abatement. Climate trends show warmer winters, slightly warmer and calmer summers, but minimal changes during the ice-free season in the southern basin (Reavie et al. 2017). These are factors that affect lake thermal conditions, internal loading, and algal seasonality and productivity.

Model limitations

With any modeling effort, we must consider its limitations, future iterations, and potential application to other lake management problems. The first key to this model's success is a nearly 60 yr record of P loading that exists for the Rainy River, which contributes 70% of the P load to LOW (Beak Consultants Ltd. 1990, Hargan et al. 2011). Although there are few lakes that have loading data with this level of historical detail (e.g., Nürnberg and LaZerte 2016), the model could be adapted to test alternative loading scenarios. We also recognize the limitations of historical monitoring data. For example, in our model we held other external P sources constant from 1850 to 2011 at 232 t P/yr (Hargan et al. 2011, HEI 2013). However, other sources include other tributary inputs, atmospheric deposition, shoreline erosion, and septic inputs, which were likely lower in presettlement times. Load monitoring of the Rainy River deserves similar scrutiny, as monitoring data from the 1950s–1970s were spotty, and we may be underestimating loads that were missed during periods of high runoff (J. Anderson, pers. comm.). Similarly, we must reconcile spotty monitoring data from the lake proper, which often recorded levels greater than 70 µg P/L in the 1960s, against low DI-TP estimates, which may be more indicative of spring TP values, during this period of peak loading (see also below). Other model components that could be refined include our model variables related to internal loading. We fix our internal loading at 2.5% of the active pool of P annually. However, if lake conditions were significantly different during the period of highest P loading (e.g., summer or winter hypolimnetic anoxia), internal loading may have historically had a greater role in nutrient dynamics. We further assume that P first entering the sediments was not buried for 10 yr, consistent with results from Lake Winnipeg sediments (Matisoff et al. 2017). Despite such model limitations and uncertainties, all combinations of variables show unequivocally that P concentrations in LOW were much higher in the past, and that the active pool of P declined over the past several decades. Most critically, we cannot create a scenario in which legacy P is a major driver of current conditions, providing a robust mechanistic argument against this hypothesis.

Management implications

Downcore profiles and model results have several important management implications for LOW and for other large shallow lake systems. First, we show that water-column concentrations of P in southern LOW declined markedly since the 1970s through nutrient abatement programs that reduced external P loading. Analysis with dynamic modeling indicates that the active pool of P was rapidly depleted from its mid-20^th century maximum via burial and outflow, and the lake has recently or should soon reach a new steady state in the absence of future stressors. The combined losses of P through outflow and burial are substantial in LOW, making the lake responsive to future reductions in external P inputs, if further loading reductions are possible. In contrast, lakes with long residence times and/or slow sedimentation rates are hampered in their ability to remove P through outflow or burial and will remain management challenges (Jeppesen et al. 2005, McCrackin et al. 2016).

Second, from a biological standpoint, we cannot say that the frequency and extent of cyanobacterial blooms is greater today than in the past in LOW. Fossil pigment records indicate that cyanobacterial blooms were also a large part of the ecology of LOW in the 1950s–1970s (Reavie et al. 2017). However, we know from fossil pigments (increase in myxoxanthophyll) that the modern cyanobacterial community is different than what was present earlier. The diatoms similarly suggest a historically unique modern scenario as communities have shifted toward more eutrophic indicators in recent decades, similar to the northern LOW “disturbed” sites studied by Rühland et al. (2010), and that diatom productivity based on biogenic silica is currently at its highest recorded levels. There is no evidence of selective downcore dissolution in the cores to suggest the upcore record is biased (Reavie et al. 2017).

It is the cause of recent algal community shifts and potential limnological shifts that must concern lake managers. Could the algal communities be responding to drivers other than P in light of the well-documented decreases in P loading and depletion of the legacy sediment P pool? Three potential drivers should be explored. Nutrient loading from the Rainy River has shifted from continuous loading to pulsed (seasonal) loading following nutrient abatement efforts that targeted sanitary and industry sources (J. Anderson, pers. comm.). Modern loadings are now highest in April–June and may have changed algal ecology where large and heavily silicified diatoms are favored in spring whereas cyanobacteria and smaller centric diatoms are favored later in the season. This response may be exacerbated by the second driver, a shift from P-limitation in the main body of LOW in the 1960s to N-limitation or co-limitation since the 1990s (Reavie et al. 2017) based on DIN:TP (Bërgstrom 2010). Although not a perfect predictor of cyanobacterial dominance (Downing et al. 2001), N-limitation has been linked to bloom formation (Ferber et al. 2004, Orihel et al. 2012).

Finally, climate changes are already evident in LOW. In its northern basins, the ice-free season has been extended by nearly 4 weeks since the 1960s (Rühland et al. 2010) with winter and summer temperatures at Kenora (Ontario) 2.3 C and 1.2 C warmer since 1900, respectively. This has resulted in increases in algal production (Paterson et al. 2017) and changes in diatom and chironomid assemblages (Rühland et al. 2008, 2010; Hyatt et al. 2011, Summers et al. 2012) that are consistent with changes in lake physical properties and water column nutrient cycling (e.g., internal loading). In contrast, the southern basin shows no discernable trend in ice-out date (MNDNR-SCO 2016). Nevertheless, climate drivers will affect the physical, chemical, and biological limnology of the lake through longer growing seasons, seasonality of external loads, and increased potential for short-term stratification. Understanding the links between these drivers, water quality, and algal ecology should be the focus of research, monitoring, and modeling on Lake of the Woods.

Acknowledgments

We thank the Minnesota Pollution Control Agency (MPCA; Contract 41642) and the Lake of the Woods Sustainability Foundation for providing funding for this project. Fieldwork was supported by the MPCA, the Sportsman's Lodge, and Captain Randy Beebe from WolfsHead Research Logistics aboard the R/V Arctic Fox. Devin Hougardy and Aaron DeRusha provided field assistance. Kathryn Hargan (PEARL, Queens University), Geoff Kramer (RESPEC), and Jesse Anderson (MPCA) provided historical loading data, Brittany Store (MPCA) provided GIS data, Matt DeWolfe (lwcb.ca) provided outflow data for the Winnipeg River. Norman Andresen (Andresen LLC, Ypsilanti, MI) analyzed diatoms in cores. Erin Mittag, Michele Natarajan, Erin Mortenson, and Alaina Fedie of the St. Croix Watershed Research Station coordinated lab analyses. Reviewers provided valuable comments that significantly improved the project.