Sediment cores

Three lake sediment cores were retrieved from Lake Superior using a Kullenberg piston corer aboard the R/V Blue Heron (Fig. 1 in text; latitude and longitude listed in Supplementary Table 1); one in 2002 (BH02-10P; Split Rock—SR), one in 2009 (BH09K-SUP09; Keweenaw—KW) and one in 2011 (BH11IR-SUP11; Isle Royale—IR). Gravity cores were also taken as trigger cores corresponding to each of the three piston cores, providing a more complete sediment core record. Eight sediment multicores, designed to capture the sediment water interface, were collected aboard the R/V Blue Heron and one aboard the R/V Lake Guardian from each basin of Lake Superior (Fig. 1 in text; latitude and longitude listed in Supplementary Table 2); one in 2003 (BH03-3), three in 2009 (BH09-2, BH09-3 and BH09-4) and four in 2010 (LG MC, IR MC, EM MC and CM MC). Core sites were identified using a Knudsen 12 kHz Hi-Res echo sounder and a CHIRP seismic reflection profiler onboard the R/V Blue Heron. The coring sites were chosen based on evidence of uninterrupted stratigraphy resulting from continuous sedimentation to characterize representative locations in Lake Superior. The LG MC core was collected aboard the R/V Lake Guardian using an Ocean Instruments model 750 box corer (30 cm × 30 cm × 90 cm), from which two 6.5 cm internal diameter cylindrical cores were sub-sampled. The rest of the multicores were collected with an Ocean Instruments model MC-400 multi-corer (9.4 cm diameter) from water depths listed in Supplementary Table 2. Initial core description and splitting were performed at the National Lacustrine Core Repository (LacCore). Subsamples were taken from each multicore at 0.5 cm intervals to 10 cm and every 1 cm thereafter to the base of the core and from each piston and accompanying gravity (trigger) core at 10 cm intervals to the base of the core. Sediment samples were subsequently frozen and freeze-dried for geochemical analyses.

Sample preparation

Homogenized dry sediment samples of known weight (12–20 mg) were placed into silver (Ag) foil capsules in designated sample trays and 1 μl of nanopure filtered H 2 O was dispensed into each sample capsule. Sample trays were then placed in a desiccator, sans desiccant, with 12 M HCl for 8 h to remove inorganic carbon from the sediment (acid fumigation). Upon removal, samples were allowed to off-gas residual HCl while being dried on a hotplate set at 60 °C; drying time was on the order of 12 h. In order to ensure dryness, samples were additionally placed in an oven set at 60 °C overnight. Dried samples were then folded in tin (Sn) foil capsules.

Elemental and isotopic analyses

TOC and total nitrogen abundances of bulk sedimentary OM were analysed for weight percent C and N concurrently with isotope ratio determinations. All concentrations are presented for acid fumigated (decarbonated) samples, eliminating variability caused by carbonate concentrations within samples. N 2 and CO 2 peak areas (Isodat v3.0) were converted to weight percent compositions using response factors generated from standards of known composition (acetanilide, caffeine, B-2153, B-2159 and urea), which were run between every ten samples. Analytical precision, based on replicate standard runs for bulk measurements, was typically better than ±0.70% on average for carbon and nitrogen weight % values and better than ±0.25‰ on average for carbon isotope values. Reproducibility between duplicate samples was better than ±0.05% for both carbon and nitrogen weight % values; reproducibility of carbon isotopic values between duplicates was better than ±0.07‰. Analyses were performed at the Large Lakes Observatory (LLO) Stable Isotope Lab using a Costech Elemental Analyzer coupled with a ThermoFinnigan DeltaplusXP stable isotope ratio monitoring mass spectrometer (EA-IRMS). Every tenth sample was run in duplicate. Carbon isotopic (δ13C) values are reported relative to Vienna Pee Dee Belemnite in conventional delta notation as per mil (‰) deviations.

Suess correction

To account for the change in the δ13C of atmospheric CO 2 from anthropogenic fossil fuel burning (the Suess effect) and consequent influence on bulk sediment δ13C values13,36 all δ13C values were corrected using the equation from ref. 36, as it encompasses the entire time period spanned by the multicores (1700–2010 AD) in this study. For lack of quantitative information on the existence of reservoir effect in Lake Superior, the correction assumed no significant time delay for particulate organic carbon sedimentation.

210Pb geochronology

Multicores were subsampled at alternating 0.5 cm intervals to 10 cm depth and alternating 1 cm intervals from 10 to 20 cm. Samples at depth were included to determine the background (supported) 210Pb activity. Unsupported 210Pb activity should be negligible at depth, where sediment age is likely to exceed 150 years (>6 half-lives). The cores were not analysed for Cs-137 because the Cs-137 peak may fall at an age younger than the time of maximum atmospheric concentration due to a lag in Cs-137 transport to its final depositional site via resuspension, transport and re-deposition in the Laurentian Great Lakes37. The 210Pb analyses for the 2003 and 2009 cores were carried out by α-spectrometry in the Department of Soil Science, University of Manitoba, under the direction of Dr Paul Wilkinson, and at Flett Research Ltd. in Winnipeg, Manitoba under the direction of Dr Robert J. Flett for the 2010 cores. The 210Pb analyses for LG MC were completed by Dr Dan Engstrom at the Science Museum of Minnesota St Croix Watershed Research Station. The age–depth relationships of the eight multicores were estimated from semi-log plots of excess 210Pb activity versus accumulated sediment mass using previously described methods38. The slopes of the straight-line segments lower in the core are proportional to sediment mass accumulation rates (MARs), which are applied to arrive at sediment age at each horizon (Supplementary Figs 1 and 2). A bioturbation zone is apparent in the uppermost part of most cores and exhibits the steepest slope38,39. The MAR for this zone is assumed to be constant, and equal to that derived from the slope of the line segment immediately below the bioturbation zone. An excursion occurs in core BH09-4 between 3 and 7 cm and is attributed to the effects of taconite ore processing beginning in 1955, with court-ordered reductions in effect by 1980; therefore, MARs for this portion are assumed to be equivalent to the slope of the segment immediately above this disturbance, but below the bioturbation zone. The 210Pb data reveal a bioturbated zone of 1.5 cm on average (Supplementary Table 2) in all cores, which is consistent with previous calculations38,39. Each zone is equivalent to 9 years in cores BH09-2, BH09-4 and BH03-3; 17, 4, 11, 13 and 18 years of sedimentation for cores BH09-3, LG MC, IR MC, EM MC and CM MC respectively. MARs in Lake Superior have not been constant over the periods of depositional history recorded in the cores—listed in Supplementary Table 2.

Paleomagnetic secular variation

Paleomagnetism is an often-used method for the age modelling of lacustrine sediments, especially where measurements of 14C are unsuccessful. Radiocarbon dates are limited in Lake Superior, as macrofossils for dating are rarely present, and sediments of gravity and piston cores present in the current study are too old for 210Pb dating. Thus, the alignment of paleomagnetic secular variation (PSV) records is the best method for dating Lake Superior sediments40. PSV records document regional variations in inclination and declination, which reflect variations in the Earth’s magnetic field with time. Ferromagnetic grains (for example, magnetite) align themselves with the local magnetic field during or shortly after deposition, upon consolidation of the sediments the motion of these grains is constrained. Any magnetization acquired by the magnetic grains long after deposition is then removed in the laboratory by alternating field demagnetization at low fields. In contrast, the natural remanent magnetization removed at higher demagnetization fields is assumed to have resulted from burial in the presence of the local magnetic field41. The PSV records from the post-glacial sediments are 10–100 year averages of the regional field, due to a 1–2 cm bioturbation zone present in the top portions of Lake Superior sediments38,39. By correlating PSV records from Superior with PSV records from regional, well-dated sites (small lakes in the region), ages can be assigned to sediment cores40. natural remanent magnetization, providing inclination and declination data for correlative purposes, was completed in 2004 at Michigan Technological University—Earth Magnetism Laboratory under the direction of Dr Suzanne Beske-Diehl for piston core BH02-10P (Split Rock—SR) and by Dr Julie Bowles at the Institute for Rock Magnetism at the University of Minnesota for the piston cores BH09K-SUP09 (Keweenaw—KW) in 2009 and BH11IR-SUP11 (Isle Royale—IR) in 2012. PSV data for the three piston cores were compared to a previously dated core, LU83-8 (ref. 42). Site-specific age–depth profiles for each of the three piston cores were completed using LU83-8 ages and associated PSV features40,42,43 (see Supplementary Fig. 3 and Supplementary Tables 3–5). For this study only features and associated ages of inclination were used in constructing the final age–depth model (regression line), due to the lack of features observable/present in the PSV declination data. Ages of inclination features apparent in each of the cores are presented as calibrated years before 1950. Correlative features for each of the three cores are listed in Supplementary Tables 3–5 and shown in Supplementary Fig. 3. The equation of the regression line (linear in SR and KW and fifth degree polynomial in IR) from the resulting plot of age versus depth from correlation with the previously dated core LU83-8 was then applied along the length of each of the cores providing the final age-depth assignments. Small amounts of Mazama ash, which is dated at 7700, cal YBP44 were identified in the Isle Royale and Keweenaw piston cores45, and included in the age models (Supplementary Fig. 4).

Age–depth models for gravity cores

Age/depth relationships for gravity cores were developed using linear extrapolation from the bottom of corresponding multicores (SR—BH09-4, KW—BH09-2 and IR—IR MC), which had previous 210Pb age assignments, as described above. Extrapolation from the 210Pb age model of the multicores to the gravity cores makes the assumption that the gravity cores begin where multicores end. Although not an ideal assumption, this is a valid approach as it provides a conservative age model, in addition to the fact that over penetration of the sediment/water interface is common when recovering both gravity and piston cores whereas multicores are designed to capture the sediment/water interface.

Diagenetic modelling

To account for the effects of diagenesis on the observed sedimentary TOC concentrations we constructed steady-state diagenetic models for each of the eight multicores, which we then compared to our measured values. Following refs 46, 47, 48, we used a reactive continuum model in which the effective bulk reactivity (k) of organic material in freshwater46 as well as marine47,48 environments decreases as a power law function of time (t) as in equation (1) where a and b are constants:

The rate of organic carbon mineralization R is defined in equation (2):

Equation (2) can be rewritten with substitution of the variable k(t) from equation (1) and subsequently integrated to obtain the fraction of organic material that remains after a given amount of time. For a≠1, separation of variables and integration of equation (3) yields equation (4), where C(t 0 ) is the concentration at the initial time t 0 .

The time t 0 corresponds to the age of OM at the sediment surface and effectively accounts for its ageing during settling through the water column and possibly for sediment resuspension46. For our purposes the initial concentration of organic carbon—C(t 0 )—is set equal to the amount of organic carbon (%TOC) measured in the uppermost sediment interval with an offset considered (listed in Supplementary Table 2). The measured TOC concentration(s) in the uppermost sediment interval may not be an accurate representation of the concentration(s) at the very interface (mathematically defined), so the offset is used to correct for potential differences. t 0 is taken to correspond to the settling time of OM in the water column (t 0 =18 years; estimated average for Lake Superior). The constants a and b were determined from the linear regression of log 10 k versus log 10 t from ref. 49 and are a=−0.985 and b=0.21.

Statistical analysis

All statistical analyses reported here were performed using SigmaPlot v. 11.0.

Generated map

Lake Superior Core Locations [map]. Scale not given. Data layers: U.S. National Park Service: World Physical Map; Esri, TomTom, U.S. Department of Commerce, U.S. Census Bureau: USA States (Generalized); mbockenhauer; Michigan Department of Technology, Management and Budget: Great Lakes Bathymetry [computer files]. University of Pittsburgh, Pittsburgh, PA: Generated by Molly D. O’Beirne, 3 April 2017. Using: ArcGIS for Desktop [GIS]. Version 10.4. Redlands, CA: Esri, 2016.

Data availability

The datasets generated during the current study are available in the PANGAEA repository (doi:10.1594/PANGAEA.874731).