Introduction

The task of restoring good water quality to anthropogenically eutrophied lakes requires solid understanding of the biogeochemical processes controlling the cycles of nutrients such as phosphorus (P) and nitrogen (N) in these systems. In many lakes, P has traditionally been considered the ultimate limiting nutrient for primary production (Schindler 1974, 1977), due to the high capacity of freshwater sediments to bind and sequester P (Blomqvist et al., 2004) and the generally poor solubility of P relative to N in terrestrial systems. Although there is ample evidence to suggest that this view is too simplistic, with many reported instances of N- or co-limitation in freshwater systems (e.g., Conley et al., 2009; Paerl et al., 2019), reduction of in-lake P concentrations remains central to numerous attempted lake restoration actions (Søndergaard et al., 2013). Hence, understanding the sources and sinks in the freshwater P cycle is important to the ongoing management of eutrophic lakes worldwide.

Multiple processes influence P cycling in lake sediments (Orihel et al., 2017). In addition to inputs of organic and inorganic P from catchment erosion, P is deposited in the form of autochthonous organic matter, which undergoes enzymatic hydrolysis in the sediment column to release dissolved phosphate (Ahlgren et al., 2005; Reitzel et al., 2007). Under oxic conditions at circum-neutral pH, phosphate may sorb to, or co-precipitate with, sedimentary Fe and Mn oxides (e.g., Slomp et al., 1996; Gunnars et al., 2002), or become incorporated into the microbial biomass as secondary organic P (Gächter et al., 1988) or polyphosphates (Ahlgren et al., 2011). Since both oxide- and polyphosphate-bound P are typically solubilized under the reducing conditions found deeper in the sediment column (Reed et al., 2011; Ahlgren et al., 2011) the phosphate thus released diffuses back towards the sediment water interface, to be incorporated into a renewed vertical cycle (Hupfer & Lewandowski, 2008). These mechanisms thus act to retain part of the accumulated P in the surface sediments, while other sediment components are buried deeper in the sediment column.

In many eutrophic lakes, external P inputs have led to significant enrichments of the surface sediments with P (Hupfer & Lewandowski, 2008, Carey & Rydin, 2011). These legacy P enrichments have been implicated in the delayed recovery from eutrophication in many systems, since part of the surface sediment enrichment may be remobilized annually to participate in water-column biogeochemical cycles. Such internal P loading may occur through reductive dissolution of the host oxide minerals during low oxygen conditions (Einsele, 1936; Mortimer, 1941), desorption of P at high pH (Koski-Vähälä et al., 2001), or physical resuspension (Evans, 1994; Niemistö et al., 2012). Due to positive feedbacks between P regeneration, enhanced primary production and anoxia (Vahtera et al., 2007), internal P loading in eutrophic systems is often significantly greater than external P inputs, complicating the management of water quality.

Understanding the timescales over which the surface sediment P enrichment may be naturally drawn down is critical to developing sensible management strategies for eutrophic lakes. Persistent internal loading of P has been suggested both empirically (Søndergaard et al., 2013) and theoretically (Katsev et al., 2006) to delay recovery from eutrophication over a timescale of decades. One of the key challenges to predicting the recovery trajectory is to estimate the rate of permanent burial of P in the sediment column. In all lakes, a certain fraction of sedimented P escapes regeneration in the surface sediments and is permanently buried and thus removed from in-lake biogeochemical cycles (e.g., Håkanson, 2003). This may include part of the originally sedimented P, as well as stable solid-phase P phases formed in the sediment column during diagenesis (Orihel et al., 2017).

A key vector of permanent P burial in lake sediments is the formation of vivianite (Fe3(PO4)2·8H2O) (e.g., Fagel et al., 2005; Rothe et al., 2016). Vivianite is a stable reduced iron (Fe) phosphate mineral formed in the sediments of lakes of various trophic states, under conditions of high availability of Fe with respect to sulfate (Rothe et al., 2015). High relative rates of dissimilatory Fe-oxide reduction with respect to sulfate reduction lead to an excess of dissolved Fe2+ in porewaters, promoting vivianite precipitation. Such conditions may be met in sedimentary microenvironments related to abandoned macrofaunal burrows (Hupfer et al., 2019), or deeper in the sediment column, for example, below the sulfate–methane transition zone (SMTZ) in brackish systems, where downwards-diffusing sulfate is exhausted by reaction with upwards-diffusing methane (Egger et al., 2015a). Due to the stability of vivianite in reduced sediments, its formation leads to permanent P burial and a drawdown of the surface sediment P enrichment over time.

By definition, anthropogenically eutrophied lakes have experienced strong perturbations in their biogeochemical cycles, leading to non-steady-state diagenesis in the sediment column (Katsev et al., 2006). The magnitude of P burial and regeneration fluxes in lakes recovering from eutrophication depend on the suite of diagenetic processes controlling vivianite formation and other P transformations in the sediment column, and how these evolve over time. Therefore, it is essential to understand the lake-specific diagenetic context in order to estimate permanent P burial in any given system. In this study, we investigate the processes controlling modern P burial and regeneration in a boreal lake recovering from a brief period of rapid eutrophication in the late twentieth century. The new information is essential to the ongoing management and restoration actions in the lake in question, and highlights key concepts applicable to the recovery of other lakes experiencing similar eutrophication trajectories.

Materials and methods

Study area and history

The Enonselkä basin forms the southernmost embayment of Lake Vesijärvi in southern Finland (Fig. 1). Enonselkä is a headwater embayment of 26.8 km2, theoretical retention time 5.6 years (Keto, 1982), fed by runoff from its catchment as well as groundwater inputs from the Salpausselkä I ice marginal formation, which lies to the south. The catchment includes forested and agricultural areas, as well as the urban area of the city of Lahti (Fig. 1). Enonselkä is generally shallow (mean depth 6.8 m) with isolated deep areas up to 30 m maximum depth. The basin is dimictic, with temperature-mediated stratification in both summer and winter. There is a net outflow of water northwards towards the adjacent Kajaanselkä basin.

Fig. 1
figure 1

Map of the Enonselkä basin of Lake Vesijärvi near Lahti, Finland. Sediment coring locations are shown in red. See Table 1 for further details

Full size image

Enonselkä was historically a clearwater system, but underwent severe eutrophication in the twentieth century due to sewage inputs from Lahti. As a consequence, severe oxygen depletion has been detected in the hypolimnion of the lake deeps since the mid-twentieth century (Keto, 1982). The construction of the Kariniemi sewage treatment plant in Lahti in 1976 led to an immediate improvement in water quality, and set the stage for the long-term recovery from eutrophication that continues to the present day (Fig. S1). However, hypolimnetic hypoxia has persisted, and the recovery has been punctuated by intervals of poorer water quality, including severe incidences of cyanobacterial blooms in the 1980s (Keto & Sammalkorpi, 1988; Keto & Tallberg, 2000). A series of restoration actions has been carried out with the goal of further improving water quality since the construction of the sewage treatment plant. These include further measures to reduce external loading, a major biomanipulation exercise in the 1990s (Horppila et al., 1998; Kairesalo et al., 1999) and most recently, aeration pumping of the hypolimnion (Niemistö et al., 2019).

Sediment coring

Sediment cores were taken from 14 sites in the Enonselkä basin (Fig. 1). A combination of short (typically 30 cm length, HTH/Kajak) and long (130–250 cm length, Livingstone) hand-held coring systems was used, and sampling was carried out over several campaigns between 2013 and 2018 (Table 1). Winter-season sampling was carried out from the ice surface, while open-water season sampling was carried out from a small motorboat. After sampling, cores were returned to laboratories in Lahti, Turku, and Helsinki for further processing. The coring sites include both shallow shelf areas, and deeps, within Enonselkä. Nine of the sites are located along three transects in the Enonsaari deep (labeled C in Fig. 1, inset). These sites were sampled within the framework of a parallel project investigating the occurrence of varved sediments in Lake Vesijärvi, and samples were made available to this study. In the present study, site W2 was investigated intensively as a case location for processes in the deeper areas of Enonselkä. Samples from other locations were included where necessary to provide supporting data.

Table 1 Sediment cores used in this study, including a summary of analyses performed. Core codes include location name (e.g., K1), as shown in Fig. 1, and HTH/Liv. to indicate the type of coring device used
Full size table

Porewater sampling

Porewaters were sampled at sites W1 and W2 (Table 1). Samples collected in February 2017 form the main body of the data presented in this study. HTH/Kajak and Livingstone core liners were pre-drilled with columns of holes for sampling by Rhizons™. The columns were taped prior to coring. At the onshore laboratory after core recovery, the tape was removed and Rhizons were inserted to extract porewater into connected 10 mL plastic syringes. Two parallel columns of holes were sampled: column I for analysis of methane (CH4) and dissolved inorganic carbon (DIC), and column II for analysis of all other porewater parameters. The syringes of the column I samples were pre-filled with 1 mL 0.1 M HNO3 to immediately acidify the sample and convert all DIC to CO2. Approximately 2 mL of porewater was extracted into each column I syringe and the exact volume was recorded for further calculations. For column II samples, up to 8 mL of porewater was extracted. The bulk sample was transferred to a centrifuge tube, from which a subsample was taken and acidified (10 µL 65% HNO3 per mL) for ICP-OES analysis. The remaining untreated bulk sample was frozen and stored for later analysis of ammonium. Bottom waters for all parameters were sampled from the overlying water in the HTH/Kajak core tube.

Porewater methane and DIC analysis

Nitrogen gas (purity 5.0) was injected into column I sample syringes to give a headspace volume of approximately 7 mL. Exact volumes for each sample were recorded. The samples were allowed to reach room temperature, then shaken to equilibrate CH4 and CO2 between headspace and dissolved phases. 5 mL of headspace gas was extracted into a dry syringe via a 3-way tap, and subsequently into a pre-evacuated 3 mL Exetainer® vial with double-wadded septum for storage (Labco Ltd., UK). Prior to analysis, helium (He) was injected into the vials to achieve a pressure of 2.0 bar as required for autoinjection into the gas chromatograph (Agilent technologies 7890B custom) equipped with a flame ionization detector (FID), electron capture detector (ECD) and thermal conductivity detector (TCD), and with He as carrier gas.

The pressurized vials were loaded into an autosampler and analyzed by gas chromatography for CH4 (FID) and CO2 (TCD). The instrument simultaneously measures N2 and O2 + Ar (TCD), allowing a 100% sum to be calculated for estimation of CH4 and CO2 concentrations in the unpressurized headspace sample (given in ppm by volume). Temporal drifts of CH4 and CO2 concentrations due to any changes in the sensitivity of detectors during long runs were accounted for by analyzing standard gas samples at 12-sample intervals as described by Penttilä et al. (2013). Headspace CH4 and CO2 concentrations were converted to total dissolved phase concentrations in the acidified porewater sample (Ctot) using Henry’s Law, as outlined for CH4 in Myllykangas et al. (2017). In Eq. 1, the first term on the right side represents the contribution to Ctot of dissolved gas released into the headspace during equilibration, while the second term represents the contribution from gas remaining in the dissolved phase:

$$C_{\text{tot}} = \frac{{X_{\text{HS}} \cdot P_{\text{atm}} \cdot V_{\text{HS}} }}{{R \cdot T \cdot V_{\text{aq}} }} + F \cdot X_{\text{HS}} \cdot P_{\text{atm}} ,$$
(1)

where XHS is the mole fraction of the gas in the headspace in ppm, Patm the pressure in the headspace in atm (= 1), VHS and VAq the syringe headspace and water volume in mL, respectively, R is the gas constant (0.08206 L atm/K/mol), T the temperature in Kelvin after equilibration (293 K), and F is the temperature-dependent equilibrium solubility coefficient of each gas in mol/L/atm, as detailed in Wiesenburg & Guinasso (1979) for CH4 and Weiss & Price (1980) for CO2. Values of Ctot were then corrected for dilution during sample acidification to give the original porewater concentration.

Porewater ICP-OES and ammonium analysis

Acidifed column II porewater subsamples were diluted with 1 M HNO3 to 10 mL volume in preparation for ICP-OES analysis. Samples were injected via a concentric nebulizer and analyzed by ICP-OES (Thermo Scientific iCAP 6000) for dissolved Fe, Mn, P and sulfur (S). Standard solution series were analyzed at regular intervals to correct for instrumental drift. Replicate analyses indicated precision of < 10%. Hydrogen sulfide, if present, is assumed to be lost during acidification, hence porewater S primarily reflects SO42− with possible minor contributions from dissolved organic S and S(0). Iron and manganese are assumed to be present in the divalent forms Fe2+ and Mn2+, respectively. Phosphorus is expected to primarily comprise ortho-phosphate (mainly as H2PO4− at in situ pH), with a minor contribution from dissolved organic P. Frozen untreated column II subsamples were thawed and analyzed for dissolved ammonium by the indophenol method (Koroleff, 1976) using an autoanalyzer (Lachat QuikChem 8000).

Sediment sampling and bulk chemical analyses

Sediment cores were sliced in preparation for chemical analysis. Sampling resolution was variable depending on the coring campaign. Samples collected in March 2013 and February 2017 from site W2 form the main body of the data presented in this study: These cores were sampled at high resolution (1, 2 or 4 cm depending on depth in core). Cores from other locations were sampled at lower resolution (see e.g., C cores in Fig. 2b). Wet sediment samples were freeze-dried or oven-dried (60°C) and water content was estimated by weight loss.

Fig. 2
figure 2

Dating procedure and estimate of linear sedimentation rates (LSR) for cores used in this study. a Dating by 137Cs of Master core C20 Liv. 2018. The 1986 peak falls within a varved interval. Varve counting confirms the Zn peak at 95 cm corresponds to 1975. b Alignment of cores using the 1975 Zn peak. c One-point age models for all cores, using the 1975 geochemical marker horizon and the sediment–water interface at the time of sampling. d Relationship between LSR and water depth for all cores. Trendline is a logarithmic least-squares fit

Full size image

Dried, homogenized sediment samples were subjected to a triple-acid digestion procedure in preparation for ICP-OES analysis. 0.1–0.2 g of sediment was extracted overnight with 5 mL HF (45%) and 5 mL of mixed HClO4 (70%)/HNO3 (65%) (volumetric ratio 3:2) at 90°C in closed Teflon vials. The acids were then evaporated at 160°C until samples displayed a gel-like consistency, and 15 mL 1 M HNO3 was added to redissolve the samples prior to analysis. ICP-OES analysis for Fe, Mn, S, P, and zinc (Zn) was then performed on the digests, in the same way as described for the porewater samples.

A subset of the dried samples from W2 Liv. 2017 was analyzed for carbon and nitrogen contents by thermal combustion elemental analysis (LECO TruSpec Micro). A subsample of 3–4 mg dried sediment was weighed into tin cups and loaded into an autosampler rosette. Inorganic carbon and nitrogen forms are considered negligible in this setting, hence no pre-extraction step was included for carbonates and measured total carbon and nitrogen contents are considered equivalent to those of the organic forms (Corg, Norg).

Dried samples from W2 Liv. 2013 were subjected to further extractions to investigate speciation of sedimentary phosphorus (P) (Table 2). Two parallel series of samples were extracted by the Bengtsson & Enell (1990) and the Hieltjes & Lijklema (1980) methods, respectively. These data were analyzed in conjunction with P data from triple-acid digested samples analyzed by ICP-OES (W2 Liv. 2017).

Table 2 Extraction protocols for the investigation of sedimentary P forms at site W2
Full size table

Sediment dating procedure

One long sediment core from the Enonsaari deep (C20 Liv. 2018) was selected on the basis of its well-preserved varved interval to be used as a reference core for sediment dating. A series of wet sediment samples from this core was analyzed for 137Cs activity using BrightSpec gamma spectrometer (3600 s counting time, Geological Survey of Finland). The single major peak in the profile was assumed to correspond to the 1986 Chernobyl nuclear accident (Ojala et al., 2017). This peak was used as the primary age marker in the core.

U-channel samples from the undisturbed surface of C20 Liv. 2018 were freeze-dried and impregnated with Araldite epoxy resin under vacuum using the method of Brauer & Casanova (2001). Varve counting was performed on the resin-embedded blocks. This process confirmed that a distinct Zn maximum observed in the ICP-OES data from this site corresponds to 1975, the last year of untreated wastewater discharge into Enonselkä. Similar Zn maxima were subsequently detected in other sediment cores and used as geochemical markers to align cores and provide a first-order estimate of linear sedimentation rate for each location. For three sites, S was used instead of Zn to identify the 1975 horizon, due to a more pronounced peak in the profile of this element (Table 1).

Sediment micro-XRF analysis

Micro-XRF analysis was performed to investigate co-location of P with other sedimentary elements in situ. U-channel samples from the uppermost 100 cm of core W2 Liv. 2017 were dehydrated with deoxygenated acetone and embedded with Spurr’s epoxy resin in preparation for Micro-XRF analysis (Jilbert et al., 2008; Jilbert & Slomp, 2013). U-channels were cut and polished to expose the internal sediment surface. The samples were then mounted in a Micro-XRF analyzer (Bruker Tornado) and a continuous series of multi-element 2D maps was produced for the depth interval 0–100 cm. Micro-XRF analysis specifications were as follows: beam diameter 25 µm, pixel resolution 25 µm, tube voltage 50 kV, anode current 599µA, measurement time 20 ms/pixel, and map width 4000 µm. Micro-XRF data were used for semi-quantitative analysis only, hence only raw count data and count ratios are reported with no further processing.

Isolation of putative vivianite crystals and SEM analysis

Seven dried sediment samples from a range of depths in core W2 Liv. 2017 were processed by a variation on the method outlined in Rothe et al. (2014) for isolation of vivianite. Briefly, freeze-dried samples were homogenized in an agate mortar, then dry-sieved to isolate the > 50 µm fraction. The > 50 µm samples were then subjected to a heavy-liquid separation procedure using lithium heteropolytungstate solution. Each sediment sample was transferred to a centrifuge tube, to which 5 mL lithium heteropolytungstate of density 2.35 g/cm3 was added. The tubes were then shaken, centrifuged, and the supernatant was removed. New solution was added and the process repeated a total of five times. The remaining heavy fraction was rinsed with deionized water and freeze-dried.

The heavy fraction samples were subsequently inspected under a reflected light microscope for evidence of vivianite. Blue-colored spherical nodules in the size range 50–150 µm (putative vivianite, as per the morphological characteristics shown in Rothe et al., 2014), when present in the samples, were isolated and mounted on double-sided tape. A selection of other components of the heavy fraction was also isolated for comparison. The tape-mounted particles were then carbon-coated and analyzed by Scanning Electron Microscopy (SEM), including Backscatter Electron (BSE) and Secondary Electron (SE) imagery, and qualitative elemental analysis by Energy Dispersive Spectroscopy (SEM–EDS).

PHREEQC determinations of vivianite saturation

The saturation index (SI) of the porewater solution with respect to vivianite (Fe3(PO4)2·8H2O) at selected depths in core W2 Liv. 2017 was calculated using PHREEQC Version 3 (Parkhurst and Appelo, 2013). A complete list of the measured and estimated input variables is given in Table S1.

Estimate of organic matter sources

A simple two end-member mixing model was applied to investigate the sources of organic matter in sediments of W2 Liv. 2017. The calculation uses only the molar N/C ratio of organic matter, and end-member values, N/CEM, based on the studies of Goñi et al. (2003) and Jilbert et al. (2018):

$$\% {\text{OC}}_{\text{phyt}} = \frac{{\left( {N/C_{\text{sample}} - N/C_{\text{EM - terr}} } \right)}}{{\left( {N/C_{\text{EM - phyt}} - N/C_{\text{EM - terr}} } \right)}} \times 100$$
(2)
$$\% {\text{OC}}_{\text{terr}} = 100 - \% {\text{OC}}_{\text{phyt}} ,$$
(3)

where N/CEM-terr = 0.04, and N/CEM-phyt = 0.13. The calculation assumes that terrestrial plant matter and phytoplankton are the only sources of organic material present in the samples, that their N/C values are spatially and temporally fixed at the end-member values, and that these values do not alter significantly during sedimentation and burial of organic matter.

Results

Age model and sedimentation rates

The signal of the 1986 Chernobyl nuclear accident is clearly visible as a peak of 2096 Bq/kg 137Cs at 69 cm in core C20 Liv. 2018 (Fig. 2a). The peak falls within a varved interval of the core, allowing the sedimentary Zn maximum at 95 cm to be constrained to the year 1975, immediately prior to the construction of the wastewater treatment plant at Kariniemi, Lahti. This geochemical marker layer is observed in all sampled cores, including the primary study core W2 Liv. 2017 (e.g., Fig. 2b), although its signal is more pronounced in S than in Zn in cores K1, K2, and K3 (Table 1). Linear sedimentation rates for the period 1975—core sampling date vary considerably between the studied sites, from 0.35 cm/yr at sites K1, K2 to 2.2 cm/yr at site C20 (Fig. 2c, d). An exponential increase in sedimentation rate is observed with increasing water depth (Fig. 3d).

Fig. 3
figure 3

Solid-phase geochemical profiles of sediments at site W2. a Bulk elemental contents by mass in core W2 Liv. 2017. Terrestrial and phytoplankton contributions to total organic carbon (Corg) were estimated according to Eqs. 2 and 3. All other elemental contents were determined by ICP-OES. b Phosphorus speciation data from core W2 Liv. 2013. Depth is adjusted to align with a. Phase 5 (Unextractable P) is set at 0.1% to account for the offset between total PTAD and total PBE (see text)

Full size image

Sediment chemical profiles at W2

The depth interval 20–60 cm in core W2 Liv. 2017, centered on the year 1975, is strongly enriched in phytoplankton-derived organic carbon, as well as P, Mn, S, and Zn (Fig. 3a). A coincident enrichment of Fe is also observed when the data are normalized to aluminum (Fe/Al). For the purposes of further discussion, we designate this interval the deep reactive layer (Fig. 3). Although the age model does not extend beyond 1975, the modern sedimentation rate at W2 of 0.93 cm/yr (Fig. 2d) suggests that the contents of these chemical components increased over several decades, before reaching a maximum in 1975. Contents are significantly reduced in the recent sediments relative to the deep reactive layer, although still elevated with respect to background values at the core base (Fig. 3a). The full 200 cm thus likely covers the entire period of urbanization of the city of Lahti, from well before 1850, through the interval of maximum eutrophication, to the modern recovery phase.

Sedimentary phosphorus speciation at W2

Of the P fractions extracted in the Hieltjes & Lijklema (1980) method, the dominant fraction throughout the upper 100 cm of W2 is 0.1 M NaOH-soluble P (Fig. 3b, Table 2). This fraction accounts for 60–70% of Total PBE, the amount of P extracted in the Bengtsson & Enell (1990) method. 1 M HCl-soluble P accounts for a further 10–15% of Total PBE, with the remainder accounted for by the whole-sample HNO3 extraction. However, Total PBE is consistently 25–30% lower than Total PTAD, the amount of P extracted by triple-acid digestion (Fig. 3a). To visualize this phase in the P speciation data, we assumed a constant background of Unextractable P (0.1% sediment by weight). With this phase included, the sum of the P species in W2 Liv. 2013 (Fig. 3b) closely matches the Total PTAD data from W2 Liv. 2017 (Fig. 3a). Total PTAD in the upper 10 cm of W2 is in the range 0.25–0.32% sediment by weight.

Microanalysis of sedimentary P at W2

Phosphorus is strongly co-located with Fe and Mn in 2D Micro-XRF maps from the deep reactive layer and shallower sediment layers at W2 Liv. 2017. Clusters of layered, mm-scale P, Fe, and Mn enrichments are clearly visible in the maps from the deep reactive layer (Fig. 5a). Heavy liquid separation of samples from this layer yielded abundant 50–150-µm-diameter blue nodules with a SEM–EDS spectrum consistent with manganous vivianite (Fig. 6a, b).

Porewater chemical profiles at W2

The deep reactive layer at W2 is characterized by convex porewater profiles of ammonium (NH4+) and DIC, up to 1.7 mmol/L and 8.0 mmol/L, respectively, at 60 cm depth (Fig. 4a). Similar convex profiles are observed in porewater Fe2+ and Mn2+ (up to 1.0 mmol/L and 0.18 mmol/L, respectively, at 60 cm), and in porewater CH4 (up to 1.8 mmol/L at 60 cm, note that this is a minimum estimate due to potential losses of methane during sampling). In the deeper sediments, the concentrations of all these parameters stabilize, although a slight decline with depth is observed in Mn2+. The 10–20 cm layer shows concave profiles with minimum values of these parameters (Fig. 4b), while higher values are observed in the uppermost layer (0–10 cm). Porewater P follows a similar pattern to NH4+ and DIC, except that P concentrations stabilize within the upper part of the deep reactive layer, reaching a maximum of 170 µmol/L, and decline towards its base (Fig. 4b). This decline continues to approximately 100 cm depth, below which concentrations rise again towards the core base. Porewater S declines from the bottom-water value of 110 µmol/L to near-zero values at 10 cm depth. The profiles of all porewater parameters at site W1 are very similar to those at W2 (Fig. S2) implying that such profiles are typical for the deeper areas of Enonselkä.

Fig. 4
figure 4

Porewater geochemical profiles from site W2, including both W2 HTH 2017 and W2 Liv. 2017. a Full (250 cm) profiles. b Upper 60 cm of the same data. Depth intervals corresponding to the deep reactive layer (20–60 cm) and zone of net oxidation (10–20 cm) are indicated. Sulfur data have been corrected for an estimated background of 20 µmol/L non-sulfate-S, hence may be considered to approximate the sulfate component only (see Fig. S6 for further details). In the zone of overlap between HTH and Livingstone cores, the lowermost four points of the HTH profile of P have been removed due to suspected analytical artefact

Full size image

Sediment chemical profiles at other sites

Phosphorus and Zn profiles similar to those at W2 are observed at all sampled sites of > 10 m water depth, implying the ubiquitous presence of the deep reactive layer in the sediment column in the deep areas of Enonselkä. The sediment depth of the 1975 maximum in P and Zn, and the absolute concentrations of these elements, increase with increasing water depth in the C sites (Fig. S3). At sites in shallower areas (i.e., 0–10 m water depth, sites labeled K in Fig. 1/Table 1), Ptot. declines from the surface downwards towards a stable background value (Fig. S4).

Discussion

Accumulation of the deep reactive layer

The profiles of many chemical components in the sediments at W2 (Fig. 3) strongly resemble the measured trajectory of electrical conductivity in Enonselkä (Fig. S1), with a rising trend towards a maximum in 1975 and a decline thereafter. The period during which the deep reactive layer was deposited was characterized by high nutrient inputs from municipal wastewater, leading to enhanced production and sedimentation of autochthonous organic matter (Keto, 1982). Urban development also led to elevated Zn concentrations in municipal and industrial wastewaters (Lahti Environment Services, personal communication), leading to elevated Zn contents in the sediments of Enonselkä. Simultaneous high water-column sulfate concentrations, related to both regional atmospheric deposition trends (Monteith et al., 2007) and local industrial sources such as pulp mills (e.g., Chen & Horan, 1998) led to high rates of sulfate reduction and sulfide mineral precipitation in the sediment column, as observed in other boreal systems (Couture et al., 2016). Finally, the enrichments of Fe, Mn, and P likely reflect enhanced redox shuttling of these elements from shallow to deep areas of the lake during intervals of widespread low oxygen conditions (Schaller & Wehrli 1996; Jilbert & Slomp, 2013; Lenz et al., 2015).

Role of the deep reactive layer in regulating modern P cycling

Although water quality in Enonselkä has improved considerably since 1975, the surface sediment enrichment of P is still far in excess of the pre-industrial background (Fig. 3), implying that full recovery of the system from the anthropogenic nutrient loading of the twentieth century remains incomplete. In the following sections, we investigate the modern diagenetic context of the sediments, in order to determine the mechanisms of P burial and regeneration before attempting to quantify their impact upon water-column P cycling in Enonselkä today.

Organic matter remineralization in the deep reactive layer

The convex profiles of ammonium (NH4+) and DIC at W2 (Fig. 4) confirm that remineralization of organic matter is presently taking place in the deep reactive layer (e.g., Burdige, 2006). At 20–60 cm depth, the layer is situated below the depth of significant bioirrigation (Kornijow & Pawlikowski, 2015) and hence microbial remineralization of organic matter is expected to proceed exclusively by anaerobic pathways. Similar convex profiles are observed in porewater Fe2+ and Mn2+, and in porewater CH4, suggesting both dissimilatory reduction of metal oxides and methanogenesis to be occurring simultaneously (Eqs. 4, 5a, 5b). Diagenetic reaction formulae used in the following sections are taken from Reed et al. (2011) and Katsev et al. (2006), with the primary redox reactions represented in the simplified forms given in Canfield (1993):

Metal oxide reduction, e.g., Fe:

$${\text{CH}}_{2} {\text{O}}_{{\left( {\text{s}} \right)}} + 7{\text{CO}}_{{2\left( {\text{aq}} \right)}} + 4{\text{Fe}}\left( {\text{OH}} \right)_{{3\left( {\text{s}} \right)}} \to 4{\text{Fe}}_{{\left( {\text{aq}} \right)}}^{2 + } + 8{\text{HCO}}_{{3\left( {\text{aq}} \right)}}^{ - } + 3{\text{H}}_{2} {\text{O}}_{{\left( {\text{l}} \right)}}$$
(4)

Methanogenesis:

$$2{\text{CH}}_{2} {\text{O}}_{{\left( {\text{s}} \right)}} + 2{\text{H}}_{2} {\text{O}}_{{\left( {\text{l}} \right)}} \to 2{\text{CO}}_{{2\left( {\text{aq}} \right)}} + 4{\text{H}}_{{2\left( {\text{aq}} \right)}}$$
(5a)
$$4{\text{H}}_{{2({\text{aq}})}} + {\text{CO}}_{{2({\text{aq}})}} \to {\text{CH}}_{{4({\text{aq}})}} + 2{\text{H}}_{2} {\text{O}}_{{({\text{l}})}}$$
(5b)

Although these reactions are typically considered to be competitive (e.g., Karvinen et al., 2015), with dissimilatory reduction of metal oxides the more energetically favorable pathway and hence observed at a shallower depth than methanogenesis (e.g., Arndt et al., 2013), vertical overlap of diagenetic zones is common in non-steady-state systems influenced by recent eutrophication (e.g., Jilbert et al., 2018; Sawicka & Bruchert, 2017). Additionally, Fe and Mn may be released into porewaters in the deep reactive layer via oxide-mediated anaerobic oxidation of methane (AOM) (Beal et al., 2009; Sivan et al., 2011; Egger et al., 2015b):

$${\text{CH}}_{{4({\text{aq}})}} + 8{\text{Fe}}({\text{OH}})_{{3({\text{s}})}} + 15{\text{H}}_{{({\text{aq}})}}^{ + } \to {\text{HCO}}_{{3({\text{aq}})}}^{ - } + 21{\text{H}}_{2} {\text{O}}_{{({\text{l}})}} + 8{\text{Fe}}_{{({\text{aq}})}}^{2 + } .$$
(6)

Diagenetic processes in the upper sediments

The porewater zonation of the upper sediments at W2 is complex and likely influenced by a combination of microbial and macrofaunal processes. The uppermost 10 cm are characterized by signals of organic matter remineralization by anaerobic pathways as described in Eqs. 4, 5a and 5b, including production of NH4+, CH4, and DIC (see zoomed profiles in Fig. 4b). The concave profile of porewater S also indicates organoclastic sulfate reduction in this layer:

$$2{\text{CH}}_{2} {\text{O}}_{{({\text{s}})}} + {\text{SO}}_{{4({\text{aq}})}}^{2 - } \to {\text{H}}_{2} {\text{S}}_{{({\text{aq}})}} + 2{\text{HCO}}_{{3({\text{aq}})}}^{ - }$$
(7)

The interval 10–20 cm is characterized by net removal of the reduced species NH4+, CH4, Fe2+, and Mn2+, as shown by concave profiles of these parameters (Fig. 4b). The most plausible explanation for these profiles is that advection of O2 by burrowing macrofauna outweighs the microbial remineralization of organic matter in this layer, leading to net consumption of reduced species in simple aerobic oxidation reactions (Eqs. 8–11). Hence we refer to this layer as the zone of net oxidation (Fig. 4b):

$${\text{NH}}_{{4({\text{aq}})}}^{ + } + 2{\text{O}}_{{2({\text{aq}})}} + 2{\text{HCO}}_{{3({\text{aq}})}}^{ - } \to {\text{NO}}_{{3({\text{aq}})}}^{ - } + 3{\text{H}}_{2} {\text{O}}_{{({\text{l}})}} + 2{\text{CO}}_{{2({\text{aq}})}}$$
(8)
$${\text{CH}}_{{4({\text{aq}})}} + 2{\text{O}}_{{2({\text{aq}})}} \to {\text{CO}}_{{2({\text{aq}})}} + 2{\text{H}}_{2} {\text{O}}_{{({\text{l}})}}$$
(9)
$$4{\text{Fe}}_{{({\text{aq}})}}^{2 + } + {\text{O}}_{{2({\text{aq}})}} + 8{\text{HCO}}_{{3({\text{aq}})}}^{ - } + {\text{H}}_{2} {\text{O}}_{{({\text{l}})}} \to 4{\text{Fe}}({\text{OH}})_{{3({\text{s}})}} + 8{\text{CO}}_{{2({\text{aq}})}}$$
(10)
$$2{\text{Mn}}_{{({\text{aq}})}}^{2 + } + {\text{O}}_{{2({\text{aq}})}} + 4{\text{HCO}}_{{3({\text{aq}})}}^{ - } \to 2{\text{MnO}}_{{2({\text{s}})}} + 2{\text{H}}_{2} {\text{O}}_{{({\text{l}})}} + 4{\text{CO}}_{{2({\text{aq}})}} .$$
(11)

The benthic macrofaunal community of Enonselkä is dominated by Chironomus and Potamothrix species (Iso-Tuisku, 2017). These organisms are capable of burrowing throughout the uppermost 20 cm of the sediment column (Hökkä, 1990) and hence advecting bottom-water oxygen into the 10–20 cm interval, as well as diluting the porewaters of this interval with bottom water. A similar influence of burrowing chironomid larvae on porewater NH4+ and Fe2+ concentrations was demonstrated by Lewandowski et al. (2007), and on porewater CH4 by Clayer et al. (2016).

Burial phases of reactive P

The observation that total PBE is consistently 25–30% lower than Total PTAD in the sediments of W2 (Fig. 3a) concurs with Hartikainen et al. (1996), who suggested that 25–30% of total P in sediments of the Enonselkä basin is bound in stable forms that cannot be extracted by most common methods. Apparently, only total sediment digestion approaches, such as our triple-acid digestion or the NaCO3 smelting method used by Hartikainen et al. (1996), are capable of extracting 100% of sedimentary P in this system. The nature of the most refractory P phase remains unclear, but it is unlikely that it participates in biogeochemical cycles due to its inaccessibility.

The large component of NaOH-soluble P (Fig. 3b) indicates that the majority of P accumulated in the upper sediments during the late twentieth century is present in reactive forms. We interpret NaOH-soluble P to represent a combination of P bound to oxides (primarily Fe- and Mn-, rather than Al-bound, as shown by Koski-Vähälä et al., 2001), labile organic P (Hartikainen et al., 1996), as well as authigenic manganous vivianite ((Fe, Mn)3(PO4)2·8H2O)). Indeed, Rothe et al. (2015) showed that vivianite is soluble in NaOH and is therefore indistinguishable from oxide-bound P in the Hieltjes & Lijklema (1980) sequential extraction (Table 2).

The Micro-XRF data allow more detailed investigation of the distribution of P in the sediments and co-location of P with other elements, and therefore can be used to augment the sequential extraction data. This approach is particularly useful for identifying heterogeneously distributed enrichments linked to authigenic or biogenic mineral formation (e.g., Jilbert & Slomp, 2013). The conspicuous clusters of P, Fe, and Mn enrichments visible in 2D Micro-XRF maps (Fig. 5a) are suggestive of an important burial phase of P at W2. To investigate how the occurrence of this phase varies with depth in the sediments, we first isolated 800 P-rich pixels from a line scan across several prominent clusters in the deep reactive layer. We then determined the P/Mn count ratio distribution of these pixels for use as a tracer of cluster occurrence throughout the upper 100 cm (Fig. 5b). P and Mn in the isolated pixels were significantly correlated to the 99% confidence level, with a least-squares regression yielding a P/Mn of 1.37 (R2 = 0.54, p = 0.88, Fig. 5b). Similarly, the median value of P/Mn in the 800 pixels was 1.33.

Fig. 5
figure 5

Micro-XRF analysis of P in resin-embedded sediments from core W2 Liv. 2017. a Example 2D map of P, Fe, and Mn in the deep reactive layer, showing co-located clusters. b Analysis of P and Mn contents (cts/s) in 800 P-rich pixels along a line scan within the deep reactive layer. c Total PTAD from ICP-OES. Depth intervals 1–7 analyzed in d are highlighted by red bars. Bars marked “x” indicate the presence of blue putative vivianite nodules isolated by heavy-liquid separation. d Analysis of P/Mn in 7 line scans at depth intervals shown in c. Each line scan contains 1000 pixels

Full size image

We subsequently plotted the frequency distribution of pixel P/Mn ratios in seven 2.5 cm scans (each representing 1000 pixels) of Micro-XRF data throughout the uppermost 100 cm of the sediments. The results show that pixel P/Mn ratios at 80–100 cm depth (scans 6 and 7 in Fig. 5c, d) are tightly distributed in the range 0–0.25, implying a relatively homogeneous distribution of P in the sediments, with a low abundance of P-Fe–Mn clusters. In the shallower layers (0–60 cm), the frequency distributions are much wider, with higher median values of P/Mn (0.4–0.7), implying an elevated contribution of P-Fe–Mn clusters to total P counts. Scan 5 (approx. 70 cm depth) is transitional between the two conditions. Hence, the appearance of the clusters roughly corresponds to the base of the deep reactive layer, and their occurrence persists upwards to the surface sediments.

The clustered co-enrichments of P, Fe, and Mn are strongly reminiscent of previous reports of authigenic manganous vivianite in freshwater and coastal marine sediments (e.g., Fagel et al., 2005, Egger et al., 2015a). The presence of vivianite is supported by two further lines of evidence. First, spherical blue nodules of 50–150 µm diameter could be isolated from the > 50 µm heavy fraction from the depths corresponding to Micro-XRF scans 1–5 in the surface sediments and deep reactive layer (Fig. 5), but could not be found in samples 6–7, which were dominated by mica and other silicates. These blue nodules bear a morphological resemblance to those presented in Rothe et al. (2014) and display SEM–EDS spectra consistent with manganous vivianite (Fig. 6a). Second, PHREEQC analysis shows that the porewater solution is supersaturated with respect to vivianite at all investigated depths (Fig. 6b). Although porewater supersaturation does not guarantee authigenic vivianite formation (Rothe et al., 2014), it is a pre-requisite for the process to occur spontaneously according to the following equilibrium reaction:

$$3{\text{Fe}}_{{({\text{aq}})}}^{2 + } + 2{\text{H}}_{2} {\text{PO}}_{{4({\text{aq}})}}^{ - } \to {\text{Fe}}_{3} ({\text{PO}}_{4} )_{{2({\text{s}})}} + 4{\text{H}}_{{({\text{aq}})}}^{ + } .$$
(12)

Hupfer et al. (2019) recently described how Fe oxide-bound P in the linings of abandoned chironomid burrows can transform in situ to vivianite in sediments of high Fe:S ratio, due to the excess of dissolved Fe over H2S during the transition to reducing conditions. This mechanism is plausible for the sediments of Enonselkä, where macrofaunal burrows exert a significant influence over the diagenetic zonation, and Fe content is strongly in excess of S (Fig. 3).

Fig. 6
figure 6

Evidence for vivianite formation in the sediments at site W2. a Spot SEM–EDS measurements on blue nodules isolated from the > 63 µm heavy-liquid separated fraction of sediments from 40 cm depth in core W2 Liv. 2017. b SEM images of a typical blue nodule, in backscatter (BSE) and secondary electron (SE) mode. c PHREEQC-derived saturation indices (SI) for vivianite in the porewaters of W2 Liv. 2017. See Table S1 for further details

Full size image

The combined evidence suggests that vivianite formation occurs in the modern sediments of Enonselkä, and therefore constitutes an important permanent burial pathway for P. Yet, part of the NaOH-soluble P fraction in the surface sediments and deep reactive layer (Fig. 3b) is very likely contributed by oxide-bound and organic P, which remains available for diagenetic reactions.

Vertical mobility of P from the deep reactive layer to surface sediments

The porewater data from W2 confirm that P is remobilized as phosphate in the upper part of the deep reactive layer (20–40 cm depth, Fig. 4). Very similar profiles were observed at site W1 (Fig. S2), implying that such mobilization from the deep reactive layer is common throughout the deep areas of the lake (> 10 m water depth). Furthermore the porewater profiles at W1 are very similar in both summer and winter (Fig. S5) indicating a continuous upwards diffusive flux from this layer throughout the year. The release of P into porewaters is expected to reflect the combined effect of direct remineralization from labile organic matter via reactions 4, 5a and 5b, as well as additional release of P associated with Fe and Mn oxides via reactions 4 and 6. This additional contribution from oxide-bound P may be significant. We calculated the theoretical porewater P profile assuming P release only from organic matter remineralization, using principles outlined in Burdige & Komada (2011):

$$\frac{{d\Delta \left[ {\text{DIC}} \right]}}{{d\Delta \left[ {{\text{H}}_{2} {\text{PO}}_{4}^{ - } } \right]}} = - r_{{{\text{C}}:{\text{P}}}} \frac{{D_{{{\text{H}}_{2} {\text{PO}}_{4}^{ - } }} }}{{D_{\text{DIC}} }},$$
(13)

where dΔ[species] = change in porewater concentration gradient over depth, rC:P = expected molar ratio of C:P in remineralized organic matter (assumed to be 106), and Dspecies = diffusion coefficients of H2PO4− and HCO3− at in situ temperature (Li & Gregory, 1974). The bottom-water P value in the calculation was fixed at the measured concentration of 15 µmol/L. This exercise shows that measured porewater P concentrations are significantly elevated with respect to the value estimated from remineralization of organic matter alone (Fig. 7), confirming that a large proportion of the P remobilized in the deep reactive layer is derived from oxide reduction.

Fig. 7
figure 7

Measured porewater P profile at site W2, Feb. 2017 (black and gray dots). In the zone of overlap between HTH and Livingstone cores, the lowermost four points of the HTH profile have been removed and the two series combined. Gray points are used to calculate the vertical diffusive flux of P from the deep reactive layer, using the measured concentration gradient of 15.02 × 10−3 µmol/cm4, assuming P to dominantly comprise H2PO4− at in situ pH of 6.5 (giving a diffusion coefficient of 4.47 × 10−6 cm2/s at in situ temperature of 4°C, Li and Gregory, 1974). White dots indicate the theoretical porewater profile expected if P would be derived exclusively from organic matter remineralization (Eq. 13). In the calculation, DIC is assumed to comprise HCO3− with a diffusion coefficient of 6.41 × 10−6 cm2/s at in situ temperature of 4°C (Li and Gregory, 1974). Colored shading indicates the deep reactive layer (grey) and zone of net oxidation (blue)

Full size image

The remobilization of P in the deep reactive layer generates a porewater concentration gradient (e.g., 15.02 × 10−3 µmol/cm4 at W2, Fig. 7), leading to a vertical upwards diffusive flux of P towards the zone of net oxidation at 10–20 cm depth. In the zone of net oxidation, in turn, the concave porewater P profile indicates net consumption of P, most likely due to co-precipitation with Fe and Mn oxides in burrows as described by Lewandowski et al. (2007). However, this immobilization of P into the solid phase does not imply that there is zero further transport towards the sediment surface. Bioturbation by Chironomus and Potamothrix is expected to mix the upper sediments, leading to biodiffusion of solid-phase constituents throughout the zone inhabited by benthic fauna (e.g., Boudreau, 1986). Part of the newly formed oxide-bound P will hence be transported vertically upwards, and taken up into repeated cycles of reduction and oxidation in the upper sediments (Hupfer & Lewandowski, 2008).

In the uppermost 10 cm, the convex porewater P profile indicates net release due to high rates of remineralization reactions. In this interval, a direct diffusive flux of P towards the sediment surface is indicated by the porewater profile (Fig. 4). When oxic conditions prevail at the sediment–water interface (in practice, for the majority of the annual cycle), the upwards diffusive flux of P in the uppermost 10 cm likely does not escape directly to the bottom water. Rather, precipitation of oxide minerals at the sediment surface traps this P in the uppermost millimeters of the sediment column, via reactions 10 and 11 (e.g., Reed et al., 2011). The actual transfer of P from the sediment to the water column in the deeper areas of the lake occurs primarily during summer hypoxic intervals, upon the reduction of oxides in the surface sediments and consequent diffusive efflux of P (Tammeorg et al. 2017). Although physical resuspension of sediments and porewaters accounts for the majority of the total flux of P from the sediments in Enonselkä (Niemistö et al., 2012), this process is primarily restricted to shallower areas (0–10 m water depth).

Basinwide estimates of P burial and regeneration from the deep reactive layer

Sedimentary P burial is the key process removing reactive P from aquatic biogeochemical cycles (e.g., Ruttenberg, 2003). The rate and distribution of P burial thus exerts a strong influence over the recovery trajectory in eutrophied aquatic systems (Jilbert et al., 2015). A typical approach to estimate permanent P burial in lake sediments under steady-state conditions is to calculate the mass accumulation rate of P at the depth at which solid-phase P contents stabilize below the active surface layer (Håkanson, 2003), often taken to be 10–25 cm (Søndergaard et al., 2003). In the deep basins of Enonselkä, this approach must be modified because of the influence of the deep reactive layer. At sites where the deep reactive layer is evident in the sediment data (in practice all sites > 10 m water depth, labeled C and W in Fig. 1/Table 1), we applied the principle that the decline in porewater P at the base of the deep reactive layer represents the ultimate lock-in depth for P in these sediments, since the diffusive gradient is reversed and the layer lies below the depth of bioturbation. Hence, solid-phase P passing through this horizon is assumed to be permanently buried. The modern burial rate of P at these sites is thus calculated by

$${\text{MAR}}_{\text{P}} = {\text{MAR}}_{{{\text{tot}}.}} \times P_{\text{tot - LID}} ,$$
(14)

where MARP = Mass accumulation rate (burial rate) of P in mg/cm2/yr, MARtot. = total mass accumulation rate of solid material in g/cm2/yr, and Ptot-LID = total PTAD content at the lock-in depth in mg/g. MARtot., in turn, was calculated from linear sedimentation rates (LSR, cm/yr) assuming volumetric porosity of 0.9 and solid-phase density of 2.65 g/cm3:

$${\text{MAR}}_{{{\text{tot}}.}} = {\text{LSR}}_{{{\text{tot}} .}} \times 0.265.$$
(15)

Burial rates of P were subsequently converted to units of mg/m2/d for comparison with earlier studies. For most sites, lock-in depth was estimated on the basis of the solid-phase profiles due to the absence of porewater data (Fig. S3). Here, the lock-in depth was assumed to be equivalent to the base of the deep reactive layer, as observed at W2.

Sites of water depth 0–10 m display more classical profiles of sedimentary P in eutrophic lakes (Fig. S4) as described in Carey & Rydin (2011). The absence of a P enrichment corresponding to the deep reactive layer at these sites implies either (1) that resuspension and focusing of sediments in these areas prevented the formation of such a layer entirely, or (2) that diagenetic processes have effectively recycled the reactive P back the sediment surface over recent decades. For these sites, we used the stable background value of Ptot. to estimate burial rates of P:

$${\text{MAR}}_{\text{P}} = {\text{MAR}}_{{{\text{tot}}.}} \times P_{{{\text{back}}.}} ,$$
(16)

where Pback. is the stable background value at the base of the core (Fig. S4).

The exponential increase in sedimentation rate observed with increasing water depth (Fig. 3d) is consistent with earlier studies reporting intensive sediment focusing into the deeps of Enonselkä (e.g., Nykänen et al., 2010; Niemistö et al., 2012; Tammeorg et al., 2018). This phenomenon strongly influences the spatial distribution of P burial (Fig. 8). To estimate basin-wide P burial rates, we first drew a logarithmic trendline through the P burial data for all sites, plotted against water depth (Fig. 8b). We then divided the basin into four water-depth classes as defined by the available bathymetric data, and used the trendline to estimate mean burial rates in each depth class. Subsequently, we multiplied these values by the areal coverage of each depth class as estimated from Fig. 1, to yield total P burial rates by depth class, as well as total basin-wide P burial. This exercise shows that the most important discrete depth class for P burial in Enonselkä is 10–20 m water depth (Fig. 8c), due to a combination of relatively high P burial rates and significant areal coverage. The total rate of P burial in Enonselkä as estimated by this method is 37,200 kg/yr, or 3.91 mg/m2/d for the entire basin on an area-normalized basis.

Fig. 8
figure 8

a Areal extent of depth classes in Enonselkä basin (as defined by map in Fig. 1). b Areal P burial rates at each coring site according to Eqs. 14 (W, C sites) and 16 (K sites). c Total P burial rates per depth class. Mean P burial rates per depth class were estimated from the logarithmic best fit line in b, and multiplied by the areal extent of depth classes shown in a

Full size image

To estimate the basin-wide upwards vertical flux of P from the deep reactive layer to the zone of net oxidation, we assumed all sediment areas below 10 m water depth (in total 5.85 km2, Figs. 1, 8a) to display similar porewater zonation to that observed at W1 and W2. Based on the porewater P profiles, the upwards flux at W1 is 0.40 mg/m2/d, and at W2, 1.34 mg/m2/d (Fig. S7). Using an average of these values (0.87 mg/m2/d) and assuming 22% of Enonselkä sediments fall in the depth range > 10 m water depth (Fig. 8a), we calculate a total upwards diffusive flux in Enonselkä of 1860 kg/yr, or 0.19 mg/m2/d for the entire basin on an area-normalized basis.

Enonselkä P budget and the relevance of the deep reactive layer for the recovery from eutrophication

We constructed a new contemporary P budget for Enonselkä basin, taking into account sediment and porewater fluxes from this study, as well as previous estimates of other components of the P cycle (Fig. 9). The budget is based on the conceptual premise that the surface sediment layer (0–20 cm) is the primary reservoir of P in Enonselkä, and that the balance of vertical fluxes in and out of this reservoir determines the trajectory and speed of trophic change in the ecosystem as a whole. We divide the surface sediment box into shallow-water areas, in which only vertical downward (i.e., burial) fluxes from the base of the surface sediment layer are permitted; and deep-water areas, in which bidirectional exchange at the base of the surface sediment layer is permitted, due to regeneration of P from the deep reactive layer. In deep areas, permanent burial occurs at the base of the deep reactive layer. Release of P from the sediments to the water column proceeds by resuspension in shallow areas, as shown by Niemistö et al. (2012), and diffusive efflux in deep areas susceptible to low oxygen conditions, as shown by Tammeorg et al. (2017). Gross sedimentation of newly produced organic P, and settling resuspended material, deposits P to the sediments throughout the lake.

Fig. 9
figure 9

Revised modern-day P budget for Enonselkä basin. All values are given to three significant figures based on available estimates, but should be considered approximate due to inherent uncertainties. Values in bold red font represent surface sediment reservoirs in kg, divided into areas of 0–10 m water depth and > 10 m water depth. Values in blue font represent fluxes in kg/yr. Water-column values are derived from Niemistö et al. (2012), corrected for the duration of the open-water season (assumed to be 7 months). Sediment values are derived from the present study, with the exception of the diffusive flux from surface sediments to water column (derived from Tammeorg et al., 2017). Sediment burial fluxes are estimated from cumulative values in Fig. 8. Sediment reservoirs are estimated assuming average thicknesses of the surface sediment and deep reactive layers as shown in the diagram, mean Ptot values of 1500 ppm (surface sediment, shallow), 2750 ppm (surface sediment, deep) and 4000 ppm (deep reactive layer), and constant porosity of 0.9

Full size image

The budget reveals that annual cycling of P between the surface sediments and water column is intense, relative to the inputs and outputs from the system as a whole. Resuspension and gross sedimentation are approximately in balance, and these fluxes cycle the equivalent of 6% of the surface sediment P pool annually. The large magnitude of these fluxes is maintained by high P contents in the surface sediments, as well as the generally shallow bathymetry of Enonselkä, which makes the system vulnerable to wind-driven mobilization of the surface sediments.

Two key further observations may be drawn from Fig. 9. Firstly, that the budget as a whole is not in steady state. There is a net annual loss of 31,150 kg P, principally through sedimentary P burial. This value is equivalent to 1.2% of the surface sediment reservoir annually, indicating that the reservoir is in steady decline. We suggest an error margin of ± 0.2% (of the surface sediment reservoir) for this value, although in practice this margin is difficult to determine without a denser network of sediment cores to better constrain the relationship in Fig. 8b. This steady loss of recyclable P concurs with the long-term trajectory towards improved water quality in Enonselkä (Fig. S1). The second observation is that the upwards flux of P from the deep reactive layer to the surface sediments is significant (1860 kg/yr), and equivalent to > 40% of the present-day value of external P loading to Enonselkä. Therefore, the deep reactive layer still plays a role in maintaining high P contents in the surface sediments today, and may be partly responsible for the slow recovery trajectory over the period 1975–present. Moreover, the present-day value for the upwards diffusive flux of 1860 kg/yr is likely to be in decline, since the porewater concentration gradient is becoming shallower as the deep reactive layer is buried. Hence, the significance of the flux may have been greater in previous decades.

Implications and conclusions

While there are examples of successful lake restoration as a consequence of reduced external P inputs, internal P loading persists in many cases (Søndergaard et al., 2013), highlighting the need to improve our understanding of P cycling below the sediment water interface in eutrophic lakes. To our knowledge, this is the first study to explicitly demonstrate the influence of a deep diagenetically active sediment layer on the P budget of a lake in recovery from eutrophication. Our results suggest that more extensive sediment and porewater profiling than hitherto attempted may be required to understand P cycling in other systems, if diagenetic processes occurring up to 60 cm deep in the sediment column could have an impact on surface sediment chemistry and P remobilization.

Although our results are of obvious use in the ongoing management of the Enonselkä basin, the concept of a deep reactive layer may be applicable in many lakes experiencing a similar eutrophication and recovery trajectory. In particular, diagenetic processes in small lakes with high sedimentation rates and historical point-source P inputs are most likely to resemble those in Enonselkä, since these conditions favor the rapid accumulation of strongly P-enriched sediment layers. Our study shows that even in Fe- and Mn-rich sedimentary settings such as Enonselkä, where P burial is expected to be relatively efficient due to vivianite formation (Katsev et al., 2006; Rothe et al., 2016), remineralization processes may continue to mobilize a fraction of sedimented P from deeply buried layers and hence to delay system recovery. In settings with lower sedimentary Fe:S ratio, this fraction may be greater and the timescale of recovery extended yet further.

The principal conclusions of the study are as follows:

  • Eutrophication has caused non-steady-state diagenesis in the sediments of the Enonselkä basin, with impacts on sedimentary P cycling

  • Due to Fe- and Mn-rich diagenetic conditions in deep areas of Enonselkä, P is rapidly buried with Fe–Mn oxides and as authigenic manganous vivianite. On the basin-scale, burial is currently in excess of all input fluxes, leading to a net drawdown of 1.2% ± 0.2% of the surface sediment P reservoir and a trajectory towards recovery

  • Sediments in the deep areas contain an organic-rich deep reactive layer, deposited at the height of eutrophication in the late twentieth century. This layer is still diagenetically active, with remineralization reactions releasing P into the porewaters

  • Upwards diffusion of P from the deep reactive layer recharges the surface sediments of deep areas of the lake with P. Surface sediments are partially mixed by bioturbating macrofauna, leading to oxidation processes and retention of upwards-diffusing P. Phosphorus is remobilized from these areas by diffusive exchange during low oxygen conditions.

  • The recharge of surface sediments with P from the deep reactive layer has likely delayed the recovery from eutrophication despite reduced external inputs. However, the influence of this flux is expected to decline as the layer is further buried over the coming decades.

  • Future management strategies for Enonselkä and similar eutrophied systems should take deep sediment P fluxes into account. Understanding the balance between permanent P burial and external P loading is critical to assessing the timescale of recovery.