Mercury records covering the past 90 000 years from lakes Prespa and Ohrid, SE Europe

Abstract

The element mercury ($\mathrm{Hg}$) is a key pollutant, and much insight has been gained by studying the present-day $\mathrm{Hg}$ cycle. However, many important processes within this cycle operate on timescales responsive to centennial- to millennial-scale environmental variability, highlighting the importance of also investigating the longer-term $\mathrm{Hg}$ records in sedimentary archives. To this end, we here explore the timing, magnitude, and expression of $\mathrm{Hg}$ signals retained in sediments over the past $\sim$ 90 $\mathrm{kyr}$ from two lakes, linked by a subterranean karst system: Lake Prespa (Greece, North Macedonia, and Albania) and Lake Ohrid (North Macedonia and Albania). Results suggest that $\mathrm{Hg}$ fluctuations are largely independent of variability in common host phases in each lake, and the recorded sedimentary $\mathrm{Hg}$ signals show distinct differences first during the Late Pleistocene (Marine Isotope Stages 2–5). The $\mathrm{Hg}$ signals in Lake Prespa sediments highlight an abrupt, short-lived peak in $\mathrm{Hg}$ accumulation coinciding with local deglaciation. In contrast, Lake Ohrid shows a broader interval with enhanced $\mathrm{Hg}$ accumulation and, superimposed, a series of low-amplitude oscillations in $\mathrm{Hg}$ concentration peaking during the Last Glacial Maximum, which may result from elevated clastic inputs. Divergent $\mathrm{Hg}$ signals are also recorded during the Early and Middle Holocene (Marine Isotope Stage 1). Here, Lake Prespa sediments show a series of large $\mathrm{Hg}$ peaks, while Lake Ohrid sediments show a progression to lower $\mathrm{Hg}$ values. Since $\sim$ 3 $\mathrm{ka}$, anthropogenic influences overwhelm local fluxes in both lakes. The lack of coherence in $\mathrm{Hg}$ accumulation between the two lakes suggests that, in the absence of an exceptional perturbation, local differences in sediment composition, lake structure, $\mathrm{Hg}$ sources, and water balance all influence the local $\mathrm{Hg}$ cycle and determine the extent to which $\mathrm{Hg}$ signals reflect local- or global-scale environmental changes.

1 Introduction

Mercury ($\mathrm{Hg}$) is a volatile metal released into the environment from both natural and anthropogenic sources and is actively cycled between surface reservoirs (e.g. atmosphere, ocean, lakes). Emissions of $\mathrm{Hg}$ by geological processes are unevenly distributed across the Earth's surface and are generally concentrated where tectonic, volcanic, and geothermal activities are most intense (Rytuba, 2003; Edwards et al., 2021; Schlüter, 2000). Geological processes have been major drivers of variability in the global $\mathrm{Hg}$ cycle throughout Earth's history (Selin, 2009), leading to the use of sedimentary $\mathrm{Hg}$ to reconstruct periods of intense volcanism (e.g. large igneous provinces; LIPs) in Earth's geological past (e.g. Grasby et al., 2019; Percival et al., 2018). In recent times, $\mathrm{Hg}$ release associated with industrialization and the extraction and combustion of fossil fuels and natural resources (metals) has overwhelmed the natural background flux (Outridge et al., 2018; Streets et al., 2019; United Nations Environment Programme, 2018).

Figure 1

A summary diagram depicting the key anthropogenic, geogenic, biospheric, cryospheric, and lacustrine processes that could generate and modify a sedimentary $\mathrm{Hg}$ signal over 10${}^{\mathrm{1}}$- to 10${}^{\mathrm{5}}$-year timescales. The processes are abbreviated as WD (wet deposition) and DD (dry deposition). Non-filled arrows depict processes acting to increase the atmospheric $\mathrm{Hg}$ burden, and colour-filled arrows depict processes acting to influence the quantity of $\mathrm{Hg}$ stored in terrestrial reservoirs. This figure is schematic (not drawn to scale) and constructed on the basis of reviews by Bishop et al. (2020), Obrist et al. (2018), and Selin (2009).

[Figure omitted. See PDF]

Existing in the atmosphere primarily in the form of gaseous elemental mercury, $\mathrm{Hg}$ has an atmospheric lifetime of up to 2 years, facilitating its deposition far from the original source (Lyman et al., 2020). Once removed from the atmosphere, $\mathrm{Hg}$ may enter vegetation and soils, where it is cycled between reservoirs by a complex series of processes, many of which occur on timescales that exceed present-day monitoring (Fig. 1) (Branfireun et al., 2020; Selin, 2009). Evasion back to the atmosphere, consumption by living organisms, or sequestration within aquatic sediments all represent ways in which $\mathrm{Hg}$ may “leave” the terrestrial environment, and aquatic sediments are known to be particularly effective sinks within the global $\mathrm{Hg}$ cycle (Bishop et al., 2020; Selin, 2009). Here, microbial processes lead to the formation of methylmercury (MeHg), which is the most bio-accumulative $\mathrm{Hg}$ species and can cause severe neurological and physiological damage to complex organisms if ingested (Driscoll et al., 2013; Wang et al., 2019).

The ecological and societal risks of environmental $\mathrm{Hg}$ contamination underscore the importance of quantifying how natural and anthropogenic processes may influence $\mathrm{Hg}$ sequestration within aquatic systems and the timescales upon which they are effective. Time-resolved sediment records sourced from marine and lacustrine basins are highly suitable for assessing these roles further back in time, as the deposited $\mathrm{Hg}$ may originate from one of several potential sources in the atmospheric (e.g. precipitation, dust), terrestrial (e.g. soils, detrital matter), aquatic, and/or lithospheric domain (Fig. 1). Thus, they can provide time-resolved records of $\mathrm{Hg}$ deposition, cycling, burial, and accumulation relative to changing environmental conditions on a local, regional, or even global scale (Cooke et al., 2020; Zaferani and Biester, 2021) and so can offer new insights into the cycling of $\mathrm{Hg}$ in the terrestrial realm.

Analysis of pre-industrial marine and lacustrine sediment records suggests that $\mathrm{Hg}$ concentration broadly reflects variability in climate (Li et al., 2020). On orbital ($>$ 10${}^{\mathrm{3}}$-year) timescales, oceanic $\mathrm{Hg}$ signals manifest as low-amplitude fluctuations corresponding to global-scale climate shifts from warm (interglacial) to colder (glacial) conditions, for example, due to changes in atmospheric composition (e.g. mineral dust loading) and circulation, biogeochemical cycling (Figueiredo et al., 2022), and/or ocean circulation (Figueiredo et al., 2020; Gelety et al., 2007; Jitaru et al., 2009; Kita et al., 2016). On centennial to millennial (10${}^{\mathrm{2}}$- to 10${}^{\mathrm{3}}$-year) timescales, lacustrine $\mathrm{Hg}$ signals correspond more closely to transient changes in hydrology, landscape dynamics, glacial ice, and permafrost extent on local-to-regional scales (Chede et al., 2022; Cordeiro et al., 2011; de Lacerda et al., 2017; Fadina et al., 2019; Li et al., 2023; Pérez-Rodríguez et al., 2018, 2015) (Fig. 1). Importantly, climate-associated $\mathrm{Hg}$ signals retained in lacustrine records integrate a range of processes, and some records show higher sedimentary $\mathrm{Hg}$ concentrations during cold, arid conditions (e.g. Li et al., 2020), while other records tend to have higher $\mathrm{Hg}$ concentrations with warm and wet climates. For example, increases in catchment-sourced detrital input have been proposed as the primary cause of $\mathrm{Hg}$ enrichment in temperate lakes (Pan et al., 2020; Schütze et al., 2018) and near-shore marine records (Fadina et al., 2019). Conversely, lakes located in glaciated regions may show dilution of $\mathrm{Hg}$ by the same inputs (Schneider et al., 2020). Local, site-specific factors are therefore likely to influence sedimentary $\mathrm{Hg}$ records. Yet, the combined effects of global and local processes complicate the study of how changes in the terrestrial $\mathrm{Hg}$ cycle may translate to measurable sedimentary signals and signals that are comparable between different regional or global archives.

Sedimentary $\mathrm{Hg}$ presence (or absence) at discrete intervals can be quantified using the total $\mathrm{Hg}$ concentration (${\mathrm{Hg}}_{\mathrm{T}}$) (Bishop et al., 2020; Kohler et al., 2022; Nasr et al., 2011). However, internal changes in bioproductivity, organic matter type and/or flux, sedimentation rate, pH, and redox conditions could all produce a distinct, local, transient sedimentary $\mathrm{Hg}$ enrichment without a meaningful change in the total amount of $\mathrm{Hg}$ present and/or mobile in the broader aquatic system. In light of these complexities, it has become common practice to examine total $\mathrm{Hg}$ concentration (${\mathrm{Hg}}_{\mathrm{T}}$) alongside $\mathrm{Hg}$ concentration divided by (normalized to) the concentration of various chemical species. Normalization is often applied when it can be shown that the abundance of a carrier (or “host”) phase directly impacts $\mathrm{Hg}$ content. Normalization (e.g. $\mathrm{Hg}$/total organic carbon, TOC; $\mathrm{Hg}$/total sulfur, TS) may, in those cases, then reveal broader changes in environmental $\mathrm{Hg}$ availability (Grasby et al., 2019; Percival et al., 2015; Shen et al., 2020; Them et al., 2019). Such an approach is particularly beneficial for studies typically spanning $>$ 10${}^{\mathrm{2}}$-year timescales, where the goal is to isolate the effects of catchment-scale depositional and/or transport processes on $\mathrm{Hg}$ signals recorded in the sediment through time.

Organic matter (hereafter represented by total organic carbon; TOC) is generally considered the dominant carrier phase of sedimentary $\mathrm{Hg}$ (Chakraborty et al., 2015; Ravichandran, 2004). For records in which TOC and $\mathrm{Hg}$ co-vary linearly, $\mathrm{Hg}$ is generally normalized to TOC (Chede et al., 2022; Figueiredo et al., 2022, 2020; Kita et al., 2016; Outridge et al., 2019). Some systems do not exhibit a relation to TOC, and $\mathrm{Hg}$ may instead be adsorbed onto (fine-grained) detrital minerals and detected by a correlation between $\mathrm{Hg}$ and mineral-dominating elements such as aluminium ($\mathrm{Al}$), titanium ($\mathrm{Ti}$), zirconium ($\mathrm{Zr}$), rubidium ($\mathrm{Rb}$), or potassium ($\mathrm{K}$) (Sanei et al., 2012; Sial et al., 2013; Them et al., 2019). In a few cases, sulfide minerals may act as important $\mathrm{Hg}$ hosts (Benoit et al., 1999; Han et al., 2008); however this is less common in freshwater lacustrine systems where sulfate reduction is often limited and only a small fraction of non-organic sulfur is buried (Ding et al., 2016; Holmer and Storkholm, 2001; Tisserand et al., 2022; Watanabe et al., 2004).

Mercury's relationship with other sedimentary components is often complex. For example, ${\mathrm{Hg}}_{\mathrm{T}}$ may also be suppressed through dilution by $\mathrm{Hg}$-poor detrital or biogenic (carbonate, silica) material, and $\mathrm{Hg}$ in many sediments is not exclusively or clearly modulated by balances between host-phase abundance and dilution. Notably, this can also occur when the host phases are always present in sufficient quantities to sequester available $\mathrm{Hg}$. In such cases, and where (single) host-phase abundance or dilution cannot be easily accounted for, the $\mathrm{Hg}$ accumulation rate (${\mathrm{Hg}}_{\mathrm{AR}}$) may provide the most optimal assessment of $\mathrm{Hg}$ availability through time as long as a robust age model is available for the archive.

Sedimentary TOC, total sulfur (TS), and detrital and biogenic mineral concentrations change in space and time, underscoring the need to assess how $\mathrm{Hg}$ covaries in relation to different host phases and other sedimentary materials. Hydrology, sedimentation regime, and geochemistry may each influence mercury host-phase availability and burial in a lacustrine system and are likely to change through time, highlighting the importance of investigating the longer-term records of $\mathrm{Hg}$ burial and accumulation in sedimentary archives.

Figure 2

(a) Map showing the location of lakes Prespa and Ohrid within southern Europe (yellow-shaded box). Volcanoes from which tephra has been identified in Co1215 (Prespa) and/or 5045-1 (Ohrid) are indicated with black triangles and numbered 1 (Vesuvius), 2 (Campi Flegrei), 3 (Ischia), 4 (Pantelleria), and 5 (Etna). Volcanoes of the South Aegean Volcanic Arc with known explosive eruptions ($>$ magnitude 4.0) between 90 and 0 $\mathrm{ka}$ are also numbered: 6 (Santorini), 7 (Nisyros), and 8 (Yali). Sites referred to in this study are also labelled as follows: (red squares) MT (Mount Tymphi), MO (Mount Olympus), MC (Mount Chelmos); (red cross) VRB (Voidomaitis River basin). (b) Aerial photo showing the coring locations of Co1215 and 5045-1 and illustrating the vegetation distributions of the area surrounding lakes Prespa and Ohrid. Mikri Prespa is labelled “MP”. Base image sourced from ©Google Earth version 9.177.0.1™. (c) Hillshade map of the Prespa–Ohrid region and bathymetric data for lakes Prespa and Ohrid presented in Jovanovska et al. (2016) and Wagner et al. (2022). The dashed grey lines denote watershed boundaries for lakes Prespa and Ohrid, adapted from Panagiotopoulos et al. (2020). Basemap sourced from ArcGIS version 10.0™ (spatial reference: 102100, 3857). The orange shadings denote mountain ranges and are labelled P/BMC (Pelister–Baba Mountain chain; the circle marks the location of Mount Pelister: 2601 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$), GMR (Galičica Mountain range), and JMR (Jablanica Mountain range; the circle marks the location of the Jablanica Mountains: 2257 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$). All mountain ranges contain evidence of the presence of glaciers and/or (peri)glacial features of the Late Pleistocene age (Hughes et al., 2022, 2023).

[Figure omitted. See PDF]

This study explores the timing, magnitude, and expression of $\mathrm{Hg}$ signals retained in the sediment records of Lake Prespa (Greece, Albania, and North Macedonia) and Lake Ohrid (North Macedonia and Albania) over the past $\sim$ 90 $\mathrm{kyr}$. The two lakes are located only $\sim$ 10 $\mathrm{km}$ apart (Fig. 2), they are hydrologically connected by karst aquifers with $\sim$ 50 % of water inflow to Lake Ohrid originating from Lake Prespa (Matzinger et al., 2006), and their sediments encode records of environmental change in southeastern Europe over the last $\sim$ 90 $\mathrm{kyr}$ (Damaschke et al., 2013; Francke et al., 2016; Leng et al., 2010; Panagiotopoulos et al., 2014; Sadori et al., 2016; Wagner et al., 2010). Comparison of their sedimentary records provides a rare opportunity to explore three important questions. First, we test how the local sedimentary environment (e.g. host-phase availability and sources) influences $\mathrm{Hg}$ burial. Second, we investigate whether $\mathrm{Hg}$ signals reflect changes in catchment hydrology, structure, and/or varying degrees of interaction between the two lake systems. Finally, we explore whether regional-scale climate variability could have measurably affected the $\mathrm{Hg}$ signals retained in the sediments.

2 Site description

2.1 Regional climate

The Mediterranean Sea and the European continent are both major influencers of present-day climate of the region surrounding lakes Prespa and Ohrid. Summer months (July to August) are hot and dry (average monthly air temperature +26 ${}^{\circ }\mathrm{C}$), while winter months (November to January) are cold, cloudy, and wet, with an average monthly air temperature of $-$1 ${}^{\circ }\mathrm{C}$ (Matzinger et al., 2006). Annual precipitation in the region averages $\sim$ 750 $\mathrm{mm}{\mathrm{yr}}^{-\mathrm{1}}$, with winter precipitation falling predominantly as snow at high elevations (Hollis and Stevenson, 1997). Present-day vegetation in the Prespa–Ohrid region comprises a mixture of endemic Balkan, central European, and Mediterranean species (Donders et al., 2021; Panagiotopoulos et al., 2014, 2020; Sadori et al., 2016).

Major shifts in the sedimentation and catchment structure of lakes Prespa and Ohrid generally correspond to the large-scale climate oscillations captured by proxy records across southern Europe throughout the last glacial–interglacial cycle ($\sim$ 100 $\mathrm{kyr}$) (e.g. Rasmussen et al., 2014; Sanchez Goñi and Harrison, 2010; Tzedakis et al., 2006). Generally higher local temperatures and moisture availability are observed during the last interglacial (pre-74 $\mathrm{ka}$), following which conditions become distinctly colder and/or drier. This resulted in the rapid recession of forest ecosystems, intense erosion of local soils and catchments, and elevated aeolian activity (e.g. Panagiotopoulos et al., 2014; Sadori et al., 2016; Francke et al., 2016). Although slightly warmer conditions were restored between $\sim$ 57 and 29 $\mathrm{ka}$, both moisture availability and temperature dropped again during the Last Glacial Maximum (LGM; $\sim$ 29–12 $\mathrm{ka}$) – favouring the growth and development of glaciers and (peri)glacial features (e.g. moraines) in the Prespa–Ohrid catchment (Ribolini et al., 2018; Gromig et al., 2018; Ruszkiczay-Rüdiger et al., 2020) but also across the Balkan Peninsula (Allard et al., 2021; Hughes and Woodward, 2017; Leontaritis et al., 2020). Lake Prespa's sediments host evidence of millennial-scale climate variability during the Last Glacial, which is tentatively correlated to Heinrich events in the North Atlantic (Wagner et al., 2010). At $\sim$ 12 $\mathrm{ka}$, the Pleistocene to Holocene transition saw the rapid propagation of warmer, wetter conditions across the region (known as Termination I), with only brief excursions from this warming trend, such as episodes of transient drying and/or cooling at 8.2 $\mathrm{ka}$ and 4.2 $\mathrm{ka}$ (Bini et al., 2019; Aufgebauer et al., 2012a). Anthropogenic influence on the Balkan landscape becomes increasingly clear from $\sim$ 2.5 $\mathrm{ka}$ onwards, mainly in the form of increased erosion regimes, forest clearance, agricultural land modification, and evidence of metallurgic practices (Panagiotopoulos et al., 2013; Cvetkoska et al., 2014; Radivojević and Roberts, 2021).

2.2 Lake Prespa

The Prespa lake system (40${}^{\circ }$54${}^{\prime }$ N, 21${}^{\circ }$02${}^{\prime }$ E) is composed of two lakes separated by an isthmus and located on the tripoint of North Macedonia, Albania, and Greece, at an altitude of 844 $\mathrm{m}$ a.s.l. The $\sim$ 1300 ${\mathrm{km}}^{\mathrm{2}}$ catchment of the Prespa lakes encompasses the Pelister Mountains to the east and the Galiçica Mountains to the southwest and west (Fig. 2). Here we focus on Megali Prespa (hereafter referred to as Lake Prespa), the larger of the two lakes, which has a surface area of 254 ${\mathrm{km}}^{\mathrm{2}}$, a maximum water depth of 48 $\mathrm{m}$, and a mean water depth of 14 $\mathrm{m}$. The total inflow into Lake Prespa averages $\sim$ 16.9 ${\mathrm{m}}^{\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$ (Matzinger et al., 2006). Water input is sourced from surface runoff (56 %), direct precipitation (35 %), and inflow from the smaller of the two lakes (Mikri Prespa; 9 %) (Matzinger et al., 2006). Lake Prespa has no surface outflow. The residence time of the lake's waters is $\sim$ 11 years (Matzinger et al., 2006), and water is predominantly lost through evaporation (52 %), underground karst channels into Lake Ohrid located 10 $\mathrm{km}$ to the west (46 %), and irrigation (2 %). The lake is currently mesotrophic with an average total phosphorus (TP) concentration of 31 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{3}}$ in the water column, basal anoxia in summer months, and generally clear waters, all signalling moderate biological productivity (Hollis and Stevenson, 1997). However, the lake likely held a more oligotrophic (low) nutrient status during the colder Late Pleistocene, during which biological productivity reduced substantially (Matzinger et al., 2006; Wagner et al., 2010).

2.3 Lake Ohrid

Lake Ohrid (41${}^{\circ }$02${}^{\prime }$ N, 20${}^{\circ }$43${}^{\prime }$ E) lies 693 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ Separated from Lake Prespa by the Galiçica Mountains, the lake straddles the boundary between North Macedonia and Albania (Fig. 2). The lake is $\sim$ 30 $\mathrm{km}$ long and 15 $\mathrm{km}$ wide, with a maximum water depth of 293 $\mathrm{m}$, water volume of 55.4 ${\mathrm{km}}^{\mathrm{3}}$, and hydraulic residence time of $\sim$ 70 years. Water input is sourced from direct precipitation (23 %), river inflow (24 %), and karst springs (53 %) fed by precipitation and water from Lake Prespa (Matzinger et al., 2006; Lacey and Jones, 2018), and this hydrological link increases the Ohrid catchment by $\sim$ 1300 ${\mathrm{km}}^{\mathrm{2}}$ to $\sim$ 2610 ${\mathrm{km}}^{\mathrm{2}}$. Evaporation (40 %) and outflow via the Crn Drim River (60 %) are the dominant ways in which water loss occurs from Lake Ohrid, and complete mixing of the lake occurs only every few years (Matzinger et al., 2006). The present-day lake shows low levels of biological productivity (oligotrophic) with an average dissolved phosphorus content of 4.5 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{3}}$, and regular mixing maintains moderately oxygenated bottom waters (Matzinger et al., 2006; Wagner et al., 2010).

3 Methods

3.1 Lake Prespa (Co1215)

Composite core Co1215 was recovered in the autumn of 2009 and summer of 2011 from the central-northern section of Lake Prespa (40${}^{\circ }$57${}^{\prime }$50${}^{\prime \prime }$ N, 20${}^{\circ }$58${}^{\prime }$41${}^{\prime \prime }$ E; Fig. 2). Sediment recovery was performed using a floating platform, with a gravity corer for surface sediments and a 3 $\mathrm{m}$ long percussion piston corer (UWITEC Co., Austria) for deeper sediments. Overlapping 3 $\mathrm{m}$ long sediment cores were cut into segments of up to 1 $\mathrm{m}$ in length for transport and storage. After splicing and correlation of core segments according to geochemical and optical information, the resulting 17.7 $\mathrm{m}$ composite core was continuously sampled at 2 $\mathrm{cm}$ resolution, yielding a total of 849 samples. It is comprised of three major lithofacies, which differ in colour, sediment structure, grain size, organic matter and carbonate content, and geochemistry. There are no lithological indications of any hiatuses or instances of non-contiguous sedimentation in core Co1215. A detailed lithostratigraphic characterization of the entire succession (90–0 $\mathrm{ka}$) is presented in Damaschke et al. (2013), along with details of the six visible tephra layers and five cryptotephra layers identified in Co1215 (Table S3 in the Supplement).

Published data for Lake Prespa (Co1215) include the following: total carbon (TC), total inorganic carbon (TIC), and total sulfur (TS) analyses (Aufgebauer et al., 2012; Damaschke et al., 2013). These data were measured at $\sim$ 2 $\mathrm{cm}$ resolution with a DIMATOC 200 (Dimatec Co., Germany) and TS using a vario MICRO cube elemental analyser (Elementar Co.) at the University of Cologne. TOC was calculated by Aufgebauer et al. (2012) as the difference between TC and TIC for the upper $\sim$ 3.2 $\mathrm{m}$ and by Damaschke et al. (2013) for the full $\sim$ 17 $\mathrm{m}$ succession. The inorganic chemistry of the sediments was determined by X-ray fluorescence (XRF) data, generated using an Itrax core scanner (Cox Ltd., Sweden) equipped with a molybdenum tube set to 30 $\mathrm{kV}$ and 30 $\mathrm{mA}$ and a $\mathrm{Si}$-drift chamber detector (Wagner et al., 2012). Core Co1215 was scanned with a resolution of 2 $\mathrm{mm}$ and a scanning time of 10 $\mathrm{s}$ per measurement. Elemental intensities were obtained for potassium ($\mathrm{K}$), titanium ($\mathrm{Ti}$), manganese ($\mathrm{Mn}$), strontium ($\mathrm{Sr}$), iron ($\mathrm{Fe}$), calcium ($\mathrm{Ca}$), and rubidium ($\mathrm{Rb}$) (Wagner et al., 2012).

Chronology

A chronology for Co1215 was previously produced by linear interpolation using volcanic ash layers, coupled with ${}^{\mathrm{14}}\mathrm{C}$ and electron spin resonance (ESR) dates obtained for bulk organic, fish, and aquatic plant remains (Aufgebauer et al., 2012). Here, we update this chronology with a Bayesian age–depth model that re-calculates previously obtained ${}^{\mathrm{14}}\mathrm{C}$ dates (Table S4 in the Supplement) with the latest (Intcal2020) radiocarbon calibration (Fig. 3) (Reimer et al., 2020). We used rBacon version 2.5.7 (Blaauw and Christen, 2011), and the new age model includes updated ${}^{\mathrm{40}}\mathrm{Ar}$/${}^{\mathrm{39}}\mathrm{Ar}$ dates of two eruptions geochemically correlated to specific tephra layers within the Prespa core (Damaschke et al., 2013): the Y-5 (39.85 $\pm$ 0.14 $\mathrm{ka}$, 2$\mathit{\sigma }$; Giaccio et al., 2017) and Y-6 (45.50 $\pm$ 1 $\mathrm{ka}$, 2$\mathit{\sigma }$; Zanchetta et al., 2018; Scaillet et al., 2013) tephra units. Every tephra layer is assumed to have been deposited instantaneously. The final model used herein presents the median of all the iterations (generally indistinguishable from the mean), and when referring to the ages of specific depths within the core we include the 95 % confidence intervals. The upper 2 $\mathrm{m}$ (Holocene) section of core Co1215 is chronologically well constrained by 10 ${}^{\mathrm{14}}\mathrm{C}$ dates and two tephra layers, with modelled age uncertainties in this section ranging from $\sim$ 5 to 580 years. Uncertainty increases with depth due to the lack of independent chronological anchors. For example, three ESR dates for a shell fragment layer ($\sim$ 14.6 $\mathrm{m}$ depth) give an average age of 73.6 $\pm$ 7.7 $\mathrm{ka}$ and form the only tie point currently available below 8.5 $\mathrm{m}$. All 27 tie points and accompanying chronological details are presented in Text S3 and Table S3 in the Supplement. Our revised model shows broad agreement with the interpolation-based chronology presented by Damaschke et al. (2013) and suggests that core Co1215 provides a continuous record of sedimentation over the past $\sim$ 90 $\mathrm{kyr}$ (Fig. S1 in the Supplement), with each 2 $\mathrm{cm}$ sample equating to $\sim$ 100 years (on average).

Figure 3

A Bayesian age–depth model for core Co1215 from Lake Prespa. The calibrated ages for the 27 tie points used in model generation are displayed by type: radiocarbon-dated bulk organic, fish, or aquatic plant remains (light grey triangles); volcanic tephra layers (black squares); and electron-spin-resonance-derived dates for a shell layer (Dreissena) located at 14.63–14.58 $\mathrm{m}$ depth (dark grey diamonds). Uncertainties for ESR dates at 1$\mathit{\sigma }$ are presented as vertical dark grey lines. The black line marks the median core age predicted by the model, which is generally indistinguishable from the predicted mean. Minimum and maximum model ages at 95 % (2$\mathit{\sigma }$) confidence are marked with orange shading. The grey bars mark the stratigraphic placement of tephra layers used as tie points, and the widths of these bars are proportional to the thickness of the tephra layers within the core. Uncertainties for radiocarbon and tephra dates are within the displayed point sizes and are presented in Table S4.

[Figure omitted. See PDF]

3.2 Lake Ohrid (core 5045-1)

The 5045-1 coring site (DEEP) is located in the central part of Lake Ohrid (41${}^{\circ }$02${}^{\prime }$57${}^{\prime \prime }$ N, 20${}^{\circ }$42${}^{\prime }$54${}^{\prime \prime }$ E) (Fig. 2). The uppermost 1.5 $\mathrm{m}$ of sediments at DEEP was recovered in 2011 using a UWITEC gravity and piston corer. Sediments below 1.5 $\mathrm{m}$ depth were recovered from six closely spaced drill holes at the site in 2013 (5045-1A to 5045-1F), with a total composite field recovery amounting to $>$ 95 % (545 $\mathrm{m}$), accounting for overlap between cores (Wagner et al., 2014b). Sediment cores were spliced to a composite record using optical and geochemical information. For sedimentological and geochemical analyses, 2 $\mathrm{cm}$ thick slices (40.7 ${\mathrm{cm}}^{\mathrm{3}}$) were removed from the core at a resolution of 16 $\mathrm{cm}$ ($\sim$ 480 years) at the University of Cologne. For this study, we analysed 217 samples between 0 and 36.27 $\mathrm{m}$ composite depth. We cannot entirely rule out that changes in sedimentation occurred between samples; however, recent seismic (Lindhorst et al., 2015), borehole logging (Ulfers et al., 2022), and sedimentological studies (Wagner et al., 2022, 2019) have suggested that sedimentation at the DEEP site has been near-continuous since $\sim$ 1.3 $\mathrm{Ma}$, with no clear evidence of any major ($>$ 1 $\mathrm{kyr}$) hiatuses. A detailed lithostratigraphic characterization of the core 5045-1 succession is presented by Francke et al. (2016). Details of the six microscopic and two visible tephra layers identified in the $\sim$ 36 $\mathrm{m}$ section analysed in this study are presented by Leicher et al. (2021) and are listed in Table S5 in the Supplement.

The $\mathrm{Hg}$ data obtained from core 5045-1 (Lake Ohrid) are presented herein alongside two previously existing datasets. The first dataset comprises TC and TIC measured using a DIMATOC 200 (TOC calculated as the difference between TC and TIC) and TS using a vario MICRO elemental analyser (Elementar Co.) at the University of Cologne – both by Francke et al. (2016). The second dataset comprises XRF data obtained using an Itrax XRF core scanner at the University of Cologne at 2.56 $\mathrm{m}$ increments, carried out on 2 $\mathrm{cm}$ thick samples and processed using QSpec 6.5 software (Cox Analytical) by Francke et al. (2016). Elemental intensities were obtained for $\mathrm{K}$, $\mathrm{Ti}$, $\mathrm{Fe}$, $\mathrm{Zr}$, and $\mathrm{Ca}$. To validate the quality of the XRF scanning data, conventional wavelength dispersive XRF (WDXRF, Philips PW 2400, Panalytical Cor.) was conducted on the 2 $\mathrm{cm}$ thick samples at 2.56 $\mathrm{m}$ resolution. Itrax data for each WDXRF sample were averaged to ensure comparability with the conventional XRF data, and $r^{\mathrm{2}}$ values were averaged to compare Itrax and WDXRF datasets (Francke et al., 2016).

Chronology

This study uses the age–depth model generated by Francke et al. (2016) and extended by Wagner et al. (2019) for the upper $\sim$ 248 $\mathrm{m}$ and $\sim$ 447 $\mathrm{m}$ of core 5045-1, respectively. Both combined tephrochronological data with orbital parameters using a Bayesian age modelling approach (Bacon 2.2). Tephra layers were used as first-order constraints. From the 11 total ${}^{\mathrm{39}}\mathrm{Ar}$/${}^{\mathrm{40}}\mathrm{Ar}$-dated tephra layers employed in Wagner et al. (2019), 7 are found in the upper $\sim$ 36 $\mathrm{m}$ section analysed in this study. The age of the eighth tie point (OH-DP-0009) is defined following geochemical correlation of this tephra layer to the 472/512 CE eruption of Somma–Vesuvius, Italy (Francke et al., 2019; Leicher et al., 2021). This chronological information was coupled with climate-sensitive proxy data (TOC and TIC) to define cross-correlation/inflection points with orbital parameters, which were included in the age–depth model as second-order constraints (Table S6 in the Supplement). Four of these points correspond to the $\sim$ 36 $\mathrm{m}$ interval analysed in this study (Wagner et al., 2019). The 95 % confidence intervals of ages for specific depths produced by the model average at $\pm$ 5.5 $\mathrm{kyr}$, with a maximum of $\pm$ 10.6 $\mathrm{kyr}$. The resulting chronology suggests that the 0.97–36.27 $\mathrm{m}$ core section analysed here covers the 1.6–89.6 $\mathrm{ka}$ time interval, with each sample possessing a resolution of $\sim$ 400 years (Francke et al., 2016; Wagner et al., 2019). A full description of the 5045-1 chronology and associated methods is presented in supplementary Text S4 in the Supplement.

3.3 Mercury measurements

Total $\mathrm{Hg}$ concentrations (${\mathrm{Hg}}_{\mathrm{T}}$) in the bulk sediments of cores 5045-1 (Ohrid) and Co1215 (Prespa) were measured using an RA-915$+$ portable mercury analyser with a PYRO-915$+$ pyrolyser, Lumex (Bin et al., 2001), at the University of Oxford. Samples were analysed for ${\mathrm{Hg}}_{\mathrm{T}}$ at a resolution of $\sim$ 2 $\mathrm{cm}$ for Co1215 (Lake Prespa) and $\sim$ 16 $\mathrm{cm}$ for 5045-1 (Lake Ohrid) (see Sects. 3.1 and 3.2). Approximately 2 ${\mathrm{cm}}^{\mathrm{3}}$ of sediment was homogenized to fine powder for TOC (Wagner et al., 2019; Francke et al., 2016; Aufgebauer et al., 2012a; Damaschke et al., 2013) and $\mathrm{Hg}$ analyses (this study). For $\mathrm{Hg}$ analysis, powdered samples were weighed into glass measuring boats, with masses ranging between 35–96 $\mathrm{mg}$ for Co1215 and between 27–78 $\mathrm{mg}$ for 5045-1. For samples particularly rich in inorganic fractions (e.g. samples coinciding with tephra layers), masses needed to be greater in order to yield a sufficiently high peak area (Lumex output) for calculation of sediment mercury concentrations. Samples were then placed into the pyrolyser (Mode 1) and heated to $\sim$ 700 ${}^{\circ }\mathrm{C}$, volatilizing any $\mathrm{Hg}$ in the sample. Spectral analysis of the gases produced the total $\mathrm{Hg}$ content of the sample. Six measures of standard material (paint-contaminated soil – NIST Standard Reference Material ^®2587) with an expected $\mathrm{Hg}$ concentration of 290 $\pm$ 9 $\mathrm{ng}{\mathrm{g}}^{-\mathrm{1}}$ (95 % confidence) were run to calibrate the instrument prior to sample analysis, and then one standard was run between every 10 lacustrine samples (calibration results are in the Supplement). Long-term observations of standard measurements with total $\mathrm{Hg}$ yield similar to the sediment samples analysed here indicate reproducibility is $\pm$ 6 % or better for $\mathrm{Hg}$ concentrations $>$ 10 $\mathrm{ng}{\mathrm{g}}^{-\mathrm{1}}$ (Frieling et al., 2023) and with $\mathrm{Hg}$ recovery close to 100 %, as expected from pyrolysis-based instrumentation for these types of rock powders (Bin et al., 2001). The details of standard runs for each core are included as a supplementary file.

Mercury accumulation

The rates of $\mathrm{Hg}$ accumulation in both cores were calculated by

1$${\mathrm{Hg}}_{\mathrm{AR}}={\mathrm{Hg}}_{\mathrm{T}}(\mathrm{DBD}\times \mathrm{SR}),$$ where ${\mathrm{Hg}}_{\mathrm{AR}}$ is the total $\mathrm{Hg}$ mass accumulation rate ($\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{kyr}}^{-\mathrm{1}}$), ${\mathrm{Hg}}_{\mathrm{T}}$ is the total mercury concentration (expressed in $\mathrm{mg}{\mathrm{g}}^{-\mathrm{1}}$), $\mathrm{DBD}$ is the dry bulk density ($\mathrm{g}{\mathrm{m}}^{-\mathrm{3}}$), and $\mathrm{SR}$ is the sedimentation rate (SR) in $\mathrm{m}{\mathrm{kyr}}^{-\mathrm{1}}$. Values for ${\mathrm{Hg}}_{\mathrm{AR}}$ are also calculated with respect to the median age estimate for each sample, meaning that uncertainties increase with depth.

Sedimentation rates for both Prespa and Ohrid were calculated by combining stratigraphic and lithological observations with the age–depth relationship ascertained for each core, respectively. For Lake Prespa, we calculate the sedimentation rate using the updated age–depth model presented in Sect. 3.1. Dry bulk density values were calculated on the basis of sedimentological data available for each core. For the Lake Ohrid dataset, DBD values were already available following the analyses of Francke et al. (2016). To acquire these values for Lake Prespa, we employed the following formula: 2$$\mathrm{DBD}=M_{\text{solid}}/V_{\text{total}},$$ where $M_{\text{solid}}$ is the mass of dry solid material ($\mathrm{g}$) measured in each sample, and $V_{\text{total}}$ is the volume of each respective sample (2 ${\mathrm{cm}}^{\mathrm{3}}$). Values for $M_{\text{solid}}$ were calculated based on recorded weight loss between wet and dry samples analysed by by Aufgebauer et al. (2012), assuming an average wet density of 1 $\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$ for wet sediments and 2.6 $\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$ (grain density) for dry sediments.

For Lake Ohrid, we utilize the sedimentation rate values calculated by Wagner et al. (2019) and dry bulk density measurements measured by Francke et al. (2016) (see these publications for full methods).

Mercury normalization

The availability of specific host phases is often assumed to exert control on the sedimentary burial of $\mathrm{Hg}$. Here, we test if the $\mathrm{Hg}$ deposited into the sediments of lakes Prespa and Ohrid may be impacted by abundance of a suite of phases. To do this, we assess both ${\mathrm{Hg}}_{\mathrm{T}}$ records relative to quantitative estimates of TOC and TS (assuming sulfides contribute to TS): both considered potential host phases of $\mathrm{Hg}$ in sedimentary successions (Chakraborty et al., 2015; Garcia-Ordiales et al., 2018; Ravichandran, 2004; Shen et al., 2020).

Detrital minerals constitute another potential host phase of $\mathrm{Hg}$ in sedimentary records. Elements such as $\mathrm{Al}$, $\mathrm{Ti}$, $\mathrm{K}$, $\mathrm{Zr}$, and $\mathrm{Rb}$ are commonly used as proxies for this purpose (Kongchum et al., 2011; Percival et al., 2018b; Shen et al., 2020). We observe a close correlation between $\mathrm{K}$ and $\mathrm{Ti}$ in Lake Prespa and quartz in Lake Ohrid (Fig. S2 in the Supplement): all proxies for fine-grained material inputs to a lake basin (Grygar et al., 2019; Warrier et al., 2016). To facilitate direct comparison of the two cores, we assess the relative abundances of (fine-grained) detrital material using XRF-based $\mathrm{K}$ counts. To account for differences in resolution between $\mathrm{Hg}$ and XRF data, $\mathrm{K}$ measurements were averaged to the thickness of each discrete $\mathrm{Hg}$ sample, and $\mathrm{K}$ values corresponding to the $\mathrm{Hg}$ sample depths were extracted.

Table 1

A comparison of the features of cores Co1215 (Lake Prespa) and 5045-1 (Lake Ohrid) relative to the Late Pleistocene (LP; 120–12 $\mathrm{ka}$), the Holocene (H; 12–0 $\mathrm{ka}$), and the marine isotope stage (MIS) stratigraphic framework defined in Lisiecki and Raymo (2005).${}^{\mathrm{a}}$ ${\mathrm{Hg}}_{\mathrm{T}}$ is given in $\mathrm{ng}{\mathrm{g}}^{-\mathrm{1}}$, and ${\mathrm{Hg}}_{\mathrm{AR}}$ is given in $\mathrm{mg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{kyr}}^{-\mathrm{1}}$.

			Depth	Mean		Sedimentology${}^{\mathrm{b}}$
			($\mathrm{m}$)	${\mathrm{Hg}}_{\mathrm{T}}$	${\mathrm{Hg}}_{\mathrm{AR}}$	Lithology	Key features
Lake Prespa	Holocene	MIS 1	2.4–0	64.6	11.9	Silt gyttja. Decreasing sand content with depth.	High lake levels. One visible and one microscopic tephra layer. High microcharcoal and green algae concentrations. High $\mathrm{TOC}/\mathrm{TN}$ ratios. High sedimentation rate.
	Late Pleistocene	MIS 2	6–2.4	41.9	12.6	2.9–2.4 $\mathrm{m}$ – high fine sand ($<$ 250 $\mathrm{\mu }\mathrm{m}$), with clayey silt and evidence of lamination.	Increasing lake level. Two cryptotephra layers. Transient nutrient pulse at 12.8–11.7 $\mathrm{ka}$. Moderate TOC and low TIC.
						6–2.9 $\mathrm{m}$ – homogenous sediment structure. Silt, distinct lamination, and siderite precipitation.	Evidence of ice-rafted debris deposition. Low productivity and lake level. High $\mathrm{K}$ and organic ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$. Low water ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$. Declining $\mathrm{C}/\mathrm{N}$ ratios. High sedimentation rate.
		MIS 3	11–6.1	32.8	7.2	6.6–6.1 $\mathrm{m}$ – massive sediment structure. Silt with distinct lamination.	Steady decrease in lake level. High oxygen index.
						11–6.6 $\mathrm{m}$ – massive sediment structure. Silt.	Increasing lake level. Four visible and three microscopic tephra layers. High $\mathrm{C}/\mathrm{N}$ ratios. Moderate TOC, very low TIC.
		MIS 4	13.9–11	33.7	9.4	Massive sediment structure. Clayey silt.	High sedimentation rate. Very low TOC. No tephra layers. Low productivity. Declining $\mathrm{C}/\mathrm{N}$ ratios. High $\mathrm{K}$ content gives evidence of ice-rafted debris deposition.
		MIS 5a–c	17.7–13.9	44.2	10.0	15.2–13.9 $\mathrm{m}$ – massive, bioturbated sediments. Clayey silt and fine sand.	Increasing lake level and high productivity. Dreissena shell layer at 14.58–14.56 $\mathrm{m}$.
						17.8–15.2 $\mathrm{m}$ – massive sediment structure. Clayey silt with fine sand.	Deep lake with moderate/low productivity. High green algae concentrations. High TOC, low TIC.
Lake Ohrid	Holocene	MIS 1	4.6–1.1	47.2	26.2	3–0 $\mathrm{m}$ – massive sediment structure. Bright colouring indicates high calcite; dark colouring indicates lower calcite.	High productivity. Four microscopic tephra layers. Low $\mathrm{K}$ concentrations. High sedimentation rate.
						4.6–3 $\mathrm{m}$ – slightly calcareous silty clay and massive sediment structure. Frequent siderite-rich layers.	Low TIC and calcite. High iron availability. Low productivity and stronger calcite dissolution. High $\mathrm{K}$ concentrations. High sedimentation rate.
	Late Pleistocene	MIS 2 MIS 3 MIS 4	11.3–4.623–11.328.8–23	69.250.650.2	45.533.429.6	Silty clay. Mottled, often massive sediment structure. Frequent siderite-rich layers. Abundant fine-fraction ($<$ 4 $\mathrm{\mu }\mathrm{m}$) sediments.	Very low TIC, TOC, and calcite, suggesting low productivity, with large inputs of fine-grained and chemically weathered siliciclastics. High iron availability. Two visible and two microscopic tephra layers. Mass-movement deposit at 7.87 $\mathrm{m}$.
		MIS 5a–c	36.3–28.8	36	20.4	35.6–28.8 $\mathrm{m}$ – silty clay with a massive sediment structure. Bright colouring indicates high calcite; dark colouring indicates lower calcite.	Low siliciclastic mineral abundance. Decreasing ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ and ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$. Strong primary productivity. Low sedimentation rate.
						36.6–35.6 $\mathrm{m}$ – silty clay. Mottled, often massive sediment structure. Frequent siderite-rich layers.	Higher carbonate ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ and ${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}$ correspond to reduced TIC and high siderite. Low sedimentation rate.

${}^{\mathrm{a}}$ MIS 5a–c, 96–71 $\mathrm{ka}$; MIS 4, 71–57 $\mathrm{ka}$; MIS 3, 57–29 $\mathrm{ka}$; MIS 2, 29–12 $\mathrm{ka}$; MIS 1, 12–0 $\mathrm{ka}$. ${}^{\mathrm{b}}$ Summarized from the following references: Lake Prespa (Aufgebauer et al., 2012; Cvetkoska et al., 2015; Damaschke et al., 2013; Leng et al., 2013; Panagiotopoulos et al., 2014; Wagner et al., 2014a); Lake Ohrid (Francke et al., 2016, 2019; Just et al., 2016; Lacey et al., 2016; Leicher et al., 2021; Wagner et al., 2019).

In line with previous studies (Shen et al., 2020), we assume that the strongest positive-sloped linear correlation with $\mathrm{Hg}$ among the analysed elements, TS, TOC, and $\mathrm{K}$, signals the most likely dominant influence on $\mathrm{Hg}$ loading in each core, which may then be interpreted as the host phase. However, it is conceivable that different host phases may dominate in different sections of the individual cores or that no single host phase clearly dominates, and so the same approach, restricted to the data within each individual marine isotope stage (MIS), is also applied (Table 1).

4 Results and discussion

Sediment cores extracted from Lake Prespa (Co1215) and Lake Ohrid (5045-1) provide a detailed, time-resolved record of $\mathrm{Hg}$ cycling between $\sim$ 90 and 0 $\mathrm{ka}$. Results are presented with direct reference to key stratigraphic intervals: the Holocene (12–0 $\mathrm{ka}$; MIS 1) and the Late Pleistocene (120–12 $\mathrm{ka}$; MISs 2–5). Widespread proxy-based evidence of warmer temperatures, forest expansion, and increased precipitation representative of interglacial climatic conditions marks the start of the Holocene epoch ($\sim$ 12 $\mathrm{ka}$) in SE Europe (Kern et al., 2022; Panagiotopoulos et al., 2014; Sadori et al., 2016; Tzedakis et al., 2006). For simplicity, we hereafter equate MIS 1 to the Holocene, allowing for a clearer distinction between glacial (Late Pleistocene) and interglacial (Holocene) climate conditions. We use these time slices, which also represent broad climate and environmental “modes”, as a framework upon which the $\mathrm{Hg}$ composition of both cores can be directly compared relative to local changes in sediment lithology and geochemistry (Table 1) and a foundation upon which local- and regional-scale environmental changes can be assessed relative to global shifts in glaciation, climate, sea level, and ocean circulation. We first consider the extent to which soft sediment processes (Sect. 4.1) and lithological features (Sect. 4.2.) may have influenced the $\mathrm{Hg}$ variability observed in Figs. 5 and 6, before adopting a catchment-scale perspective in Sect. 4.3 to explore the role of diverse environmental processes in $\mathrm{Hg}$ cycling through these two systems.

Figure 4

A comparison of host-phase relationships between lakes Prespa and Ohrid. The points are coded relative to the stratigraphic period: the Holocene (12–0 $\mathrm{ka}$, transparent circles) and the Late Pleistocene (90–12 $\mathrm{ka}$, filled symbols). We compare ${\mathrm{Hg}}_{\mathrm{T}}$ records for both lakes relative to total organic carbon (TOC), sulfide (estimated by total sulfur; TS), and detrital minerals (estimated by potassium ($\mathrm{K}$) concentrations) – note that aluminium (Al) data are more commonly used as an indicator of detrital mineral abundance, but these are currently unavailable for 5045-1.

[Figure omitted. See PDF]

4.1 Host-phase controls

The availability and abundance of specific host phases are often assumed to control sedimentary $\mathrm{Hg}$ accumulation and burial (Outridge et al., 2007). Both Lake Prespa and Lake Ohrid show evidence of complex relationships among ${\mathrm{Hg}}_{\mathrm{T}}$, TOC, TS, and $\mathrm{K}$ concentrations through time (Fig. 4). However, the trends displayed in Fig. 4 also suggest the following: (1) the strength of the relationships among $\mathrm{Hg}$, TOC, TS, and detrital minerals ($\mathrm{K}$) is distinctly different between the two lakes, and (2) the ${\mathrm{Hg}}_{\mathrm{T}}$ signals preserved in Lake Prespa and Lake Ohrid cannot be fully and individually explained by variability in abundance of these potential host phases individually.

Figure 5

Total $\mathrm{Hg}$ (${\mathrm{Hg}}_{\mathrm{T}}$) and total $\mathrm{Hg}$ accumulation rate (${\mathrm{Hg}}_{\mathrm{AR}}$) for core Co1215 from Lake Prespa, presented as a function of depth and time. Displayed are key sediment lithofacies, as well as visible (grey shading) and microscopic (orange shading) tephra layers. We include records of ${\mathrm{Hg}}_{\mathrm{T}}$ (this study) normalized to records of total organic carbon (TOC) (Damaschke et al., 2013), total sulfur (TS) (Aufgebauer et al., 2012), and detrital mineral abundance (estimated by potassium; $\mathrm{K}$) (Panagiotopoulos et al., 2014), with filled shading marking the original datasets. A distinct lake lowstand based on seismic profiles and sedimentological data is marked at 14.63–14.58 $\mathrm{m}$ depth (red shading) (Wagner et al., 2014a). The purple arrow marks sections where artificially high $\mathrm{Hg}/\mathrm{TOC}$ values are generated by a sharp drop to near-0 TOC ($<$ 0.06 $\mathrm{wt}\mathrm{\%}$), coinciding with deposition of the Y-5 (17.1 $\mathrm{m}$) tephra unit – an effect expected as background sedimentation is interrupted by volcanic ash deposition. The white boxes mark the marine isotope stages defined by Lisiecki and Raymo (2005), and stratigraphic periods are labelled in black/white.

[Figure omitted. See PDF]

Core Co1215 from Lake Prespa shows a moderate correlation between ${\mathrm{Hg}}_{\mathrm{T}}$ and TOC during the Holocene and Late Pleistocene (all data in Fig. 4; Table 1). This correlation is most significant during the Holocene (MIS 1), where distinct enrichments in ${\mathrm{Hg}}_{\mathrm{T}}$ occur in conjunction with a similarly sharp increase in TOC and low variability in $\mathrm{Hg}/\mathrm{TOC}$ values (Fig. 5). However, it is more inconsistent during the Late Pleistocene (MIS 2–5). For example, the highest ${\mathrm{Hg}}_{\mathrm{T}}$ values are measured in the relatively TOC-lean sediments of MIS 2 (Figs. 4 and 5), and a plateau also appears when higher TOC concentrations are reached during MIS 5, whereby ${\mathrm{Hg}}_{\mathrm{T}}$ no longer increases in step with TOC (Figs. 4 and S2). The correlations observed are not strong enough to conclude that TOC availability can fully explain the $\mathrm{Hg}$ signals observed in Lake Prespa throughout the $\sim$ 90 $\mathrm{kyr}$ succession.

Correlations among ${\mathrm{Hg}}_{\mathrm{T}}$, detrital mineral, and/or TS availability are also largely absent, suggesting that the complex Hg–TOC relationship is not a function of time-varying sulfides and detrital mineral availability. Large peaks in $\mathrm{Hg}/\mathrm{K}$ are visible during the Holocene (Fig. 5), but these are not reflected in ${\mathrm{Hg}}_{\mathrm{AR}}$ and are therefore an artefact of considerably lower $\mathrm{K}$ concentrations within this section of the core rather than indicators of changes in lake $\mathrm{Hg}$ levels. The highest positive $r^{\mathrm{2}}$ value between ${\mathrm{Hg}}_{\mathrm{T}}$ and TS is observed during the Holocene (MIS 1: $r^{\mathrm{2}}$ $=$ 0.25) (Fig. 4), implying that $>$ 75 % of variance in the dataset cannot be explained with sulfide availability during this time period. Correlations for other periods are even weaker, and some periods appear to show distinct patterns of $\mathrm{Hg}$ and potential host-phase behaviour (Fig. 4).

One possibility is that $\mathrm{Hg}$ signals reflect changes in the dominant sources of organic and detrital materials deposited into the lake. For example, combined isotopic and sedimentological data record episodes of stronger algal blooms during MIS 1 and 5 (Leng et al., 2013), supported by coeval abundance of freshwater diatom genera such as Cyclotella and Aulacoseira (Cvetkoska et al., 2015). All correspond to elevated ${\mathrm{Hg}}_{\mathrm{T}}$ and so could imply more effective $\mathrm{Hg}$ burial by autochthonous organic material compared to allochthonous organic material (Leng et al., 2013; Damaschke et al., 2013). However, in the presence of abundant binding ligands such as for the Lake Prespa record, maximum $\mathrm{Hg}$ burial is limited principally by supply regardless of productivity, and so changing $\mathrm{Hg}$ signals in Lake Prespa more likely reflect changes in environmental $\mathrm{Hg}$ availability, resulting from externally driven oscillations in $\mathrm{Hg}$ emission and/or exchange between (local) surface reservoirs such as forests, water courses, and soils (Bishop et al., 2020; Obrist et al., 2018). This interpretation is supported by the lack of a close statistical correspondence among $\mathrm{Hg}$, organic matter, and sulfur or detrital mineral content, source, and composition (Fig. 4), which suggests that $\mathrm{Hg}$ burial efficiency is only weakly associated with host-phase availability in this system.

Figure 6

Total $\mathrm{Hg}$ (${\mathrm{Hg}}_{\mathrm{T}}$) and total $\mathrm{Hg}$ accumulation rate (${\mathrm{Hg}}_{\mathrm{AR}}$) for core 5045-1 from Lake Ohrid, presented as a function of depth and time. Displayed are key sediment lithofacies, as well as visible (grey shading) and microscopic (orange shading) tephra layers. We include records of ${\mathrm{Hg}}_{\mathrm{T}}$ (this study) normalized to records of total organic carbon (TOC) (Francke et al., 2016), total sulfur (TS) (Francke et al., 2016), and detrital mineral abundance (estimated by potassium; $\mathrm{K}$) (Francke et al., 2016; Wagner et al., 2019), with filled shading marking the original datasets. A mass-movement deposit (MMD) is marked at 7.87 $\mathrm{m}$ depth (brown shading) (Francke et al., 2016). The purple arrows mark sections where artificially high $\mathrm{Hg}/\mathrm{TOC}$ values are generated by a sharp drop to near-0 TOC ($<$ 0.06 $\mathrm{wt}\mathrm{\%}$), coinciding with deposition of the Y-5 (17.1 $\mathrm{m}$) and Mercato (11.5 $\mathrm{m}$) tephra layers – an expected effect as background sedimentation is interrupted by volcanic ash deposition. The white boxes mark the marine isotope stages as defined by Lisiecki and Raymo (2005), and stratigraphic periods are labelled in black/white.

[Figure omitted. See PDF]

Core 5045-1 from Lake Ohrid shows elevated ${\mathrm{Hg}}_{\mathrm{T}}$ during the Late Pleistocene compared to the Holocene (Fig. 6; Table 1). Peaks in ${\mathrm{Hg}}_{\mathrm{T}}$ most consistently correspond to increases in $\mathrm{K}$ (detrital mineral) intensities, reflected by a broadly positive relationship between ${\mathrm{Hg}}_{\mathrm{T}}$ and K throughout the succession (Figs. 4 and S3 in the Supplement). However, this relationship is only described by $r^{\mathrm{2}}$ values $<$ 0.5, and the strength of this correlation varies across the span of the record, weakening during the Holocene (Fig. 4).

Variable $\mathrm{Hg}$ values in the Ohrid record appear to be less influenced by organic matter and/or sulfide availability. Fluctuations in TOC/TS values suggest that some sulfide formation may have occurred during the Late Pleistocene (MISs 2–5) (Wagner et al., 2009; Francke et al., 2016). However, even in these phases, TS remains low and correlations between ${\mathrm{Hg}}_{\mathrm{T}}$ and TS are generally negative or weak ($r^{\mathrm{2}}$ $<$ 0.2; Fig. 4) so that $\mathrm{Hg}$ signals do not change in magnitude or expression even when TS variability is accounted for (Fig. 6), potentially due to the oligotrophic state of Lake Ohrid favouring burial of sulfide-depleted sediments (Francke et al., 2016; Vogel et al., 2010). More remarkably, the relationship between ${\mathrm{Hg}}_{\mathrm{T}}$ and organic matter in Lake Ohrid also shows an inverse correlation (Fig. 4). These trends may be explained by a scenario where the $\mathrm{Hg}$ flux to Ohrid from direct deposition and/or the surrounding catchment is typically the limiting factor, rather than availability of potential host phases.

Lake Ohrid and Lake Prespa show distinct differences in the strength of their $\mathrm{Hg}$–host-phase relationships. In Lake Prespa, $\mathrm{Hg}$ broadly covaries with organic matter (TOC), whereas in Lake Ohrid correlations are observed between $\mathrm{Hg}$ and detrital minerals ($\mathrm{K}$). Nonetheless, only a relatively small proportion of $\mathrm{Hg}$ variability can be explained by host-phase availability in each record. This suggests that while host-phase availability may, at times, exert an influence on the $\mathrm{Hg}$ signals recorded in these lakes, the catchment-controlled changes in $\mathrm{Hg}$ fluxes are typically the more dominant effect on $\mathrm{Hg}$ in these sediment records. In the absence of a pronounced host-phase influence, retention of a measurable $\mathrm{Hg}$ signal requires that the net influx of $\mathrm{Hg}$ into the lake (e.g. surface runoff, wet/dry deposition) exceeds the amount leaving the system due to processes such as runoff or evasion. Therefore, we surmise that the ${\mathrm{Hg}}_{\mathrm{T}}$ and ${\mathrm{Hg}}_{\mathrm{AR}}$ signals recorded in Lake Prespa and Lake Ohrid are records of net $\mathrm{Hg}$ input into the two lakes rather than the efficiency of sedimentary drawdown.

4.2 Tephra layers

As volcanic eruptions are among the most significant natural $\mathrm{Hg}$ sources, we assess whether the previously recognized tephra deposition events in Lake Prespa correspond to changes in $\mathrm{Hg}$ deposition. Overall, we find that individual tephra horizons and surrounding sediments do not consistently correspond to measurable peaks in ${\mathrm{Hg}}_{\mathrm{T}}$ or ${\mathrm{Hg}}_{\mathrm{AR}}$ in Lake Prespa (Fig. 5). Only 2 of the 11 preserved ash layers coincide with elevated ${\mathrm{Hg}}_{\mathrm{T}}$: Mercato (8.54 $\pm$ 0.09 $\mathrm{ka}$; Somma-Vesuvius) and LN1 (14.75 $\pm$ 0.52 $\mathrm{ka}$; Campi Flegrei). These two units are not associated with disproportionately large tephra volumes, and neither coincides with evidence of transient changes in authigenic carbonate precipitation or sediment diagenesis that may impact sedimentary $\mathrm{Hg}$. This implies that $\mathrm{Hg}$ concentrations in Lake Prespa cannot, in general, be unequivocally linked to short-lived ($<$ 1 year) individual eruption events between $\sim$ 90 and 0 $\mathrm{ka}$ (Fig. S5 in the Supplement).

Discrete ash fall events (recorded by tephra/cryptotephra) do not consistently correspond to measurable peaks in ${\mathrm{Hg}}_{\mathrm{T}}$ or ${\mathrm{Hg}}_{\mathrm{AR}}$ in the slightly lower-resolution ($\sim$ 400 years per sample) Lake Ohrid record (Fig. S5). Considering this lack of correspondence of $\mathrm{Hg}$ with ash layers, in conjunction with the Lake Prespa data too, suggests that (a) surface $\mathrm{Hg}$ loading was not appreciably increased with most large eruption events over the past 90 $\mathrm{kyr}$ in the Balkans and/or (b) sampling resolution may need to be significantly higher and/or focused on less-bioturbated records to identify single, short-lived volcanogenic perturbations of the scale and type occurring during the period recorded in the Ohrid (and Prespa) sedimentary successions.

Figure 7

Total mercury (${\mathrm{Hg}}_{\mathrm{T}}$) and mercury accumulation rate (${\mathrm{Hg}}_{\mathrm{AR}}$) records for Lake Prespa and Lake Ohrid generated by this study and proxy datasets generated by prior studies. For Lake Prespa, these include arboreal pollen (AP) concentrations (Panagiotopoulos et al., 2014), microcharcoal (Panagiotopoulos, 2013), potassium ($\mathrm{K}$) (Aufgebauer et al., 2012; Wagner et al., 2010), and pollen assemblage zones (PAZs) (Panagiotopoulos et al., 2014). For Lake Ohrid, these include AP concentrations (Sadori et al., 2016); potassium ($\mathrm{K}$) (Wagner et al., 2019; Francke et al., 2016); 1000-year average surface–air temperature (SAT; ${}^{\circ }\mathrm{C}$) and annual mean precipitation (millimetres), both simulated by the LOVECLIM Earth system model (Goosse et al., 2010) for the Prespa–Ohrid region (Wagner et al., 2019); and pollen assemblage zones (PAZs) (Sadori et al., 2016). Pollen assemblage zones, defined by Panagiotopoulos et al. (2014) (Lake Prespa) and Sadori et al. (2016) (Lake Ohrid), are presented as green bars, shaded relative to the tree population density (the darker the colour, the higher the density). We include a chronology of glacial processes based on radiometric dating of glacial landforms in the following locations: the Voidomatis River basin (purple) (Lewin et al., 1991; Woodward et al., 2008), the Pindus Mountains (lilac) (Allard et al., 2021, 2020; Styllas et al., 2018; Hughes et al., 2006; Pope et al., 2017), and the Dinaric Alps (blue) (Gromig et al., 2018; Ribolini et al., 2018; Ruszkiczay-Rüdiger et al., 2020). The white boxes mark the marine isotope stages (MISs) as defined by Lisiecki and Raymo (2005), and stratigraphic periods are labelled in black/white. Vertical grey shading denotes the timing of the largest changes in glacier extent and volume.

[Figure omitted. See PDF]

4.3 Variability through time

4.3.1

Late Pleistocene (90–35 $\mathrm{ka}$; MIS 5 to MIS 3)

The Lake Prespa and Lake Ohrid sediment cores show similarly muted variability in ${\mathrm{Hg}}_{\mathrm{T}}$ and ${\mathrm{Hg}}_{\mathrm{AR}}$ values between $\sim$ 90 and 35 $\mathrm{ka}$ (broadly MISs 5a–c, 4, and 3), alluding to relatively stable $\mathrm{Hg}$ inputs (Fig. 7; Table 1). High organic and low clastic material concentrations point to warmer climate conditions during this interval, in which both catchments experienced an increase in moisture availability, pronounced forest expansion, and plant diversification – collectively acting to stabilize hillslopes and reduce deep soil erosion (Francke et al., 2019; Panagiotopoulos et al., 2014; Sadori et al., 2016). One possibility is that $\mathrm{Hg}$ sequestration during this interval was controlled by consistent rates of algal scavenging (Biester et al., 2018; Outridge et al., 2007, 2019; Stern et al., 2009). Elevated TOC (Fig. 5), hydrogen index, $\mathrm{TOC}/\mathrm{TN}$, and biogenic carbonate concentrations between $\sim$ 90 and 71 $\mathrm{ka}$ in both Lake Prespa and Lake Ohrid signal nutrient upwelling and increased allochthonous inputs, in conjunction with elevated primary productivity. For example, Lake Prespa records green algae accumulation (Cvetkoska et al., 2016, 2015; Leng et al., 2013; Panagiotopoulos et al., 2014), and sediments rich in biogenic silica (${\mathrm{bSiO}}_{\mathrm{2}}$) are also evident in Lake Ohrid (Francke et al., 2016). Slow changes in lake geochemistry associated with these biological processes are consistent with a steady ${\mathrm{Hg}}_{\mathrm{AR}}$ in both Lake Prespa and Lake Ohrid during this time and the absence of any especially pronounced changes in ${\mathrm{Hg}}_{\mathrm{T}}$. This could suggest that, for a relatively prolonged period ($\sim$ 96–35 $\mathrm{ka}$), $\mathrm{Hg}$ flux to the two lakes did not change with a magnitude sufficient to cause measurable sedimentary changes, and processes capable of amplifying differences in sedimentary $\mathrm{Hg}$ between Ohrid and Prespa were not particularly influential.

MIS 3 marks the start of slow divergence between the $\mathrm{Hg}$ records of Lake Prespa and Lake Ohrid. During MIS 3, proxy records suggest that conditions in the Prespa–Ohrid region were milder than in MIS 4 but cooler and drier than in MIS 5 (Fig. 7) (Panagiotopoulos et al., 2014; Sadori et al., 2016; Wagner et al., 2019). Divergent $\mathrm{Hg}$ signals could be linked to two climate-driven processes. First, a reduction in primary productivity in Lake Prespa signalled by decreasing TOC, hydrogen index, and endogenic carbonate compared to values observed during MIS 5 (Aufgebauer et al., 2012; Cvetkoska et al., 2016; Leng et al., 2013). Second is an increase in detrital material flux to both lakes (signalled by elevated K count; Fig. 7) due to recession of the surrounding forests and subsequently elevated rates of catchment erosion (Damaschke et al., 2013; Francke et al., 2019; Panagiotopoulos et al., 2014; Sadori et al., 2016). This environmental shift is more likely to favour enhanced $\mathrm{Hg}$ mobility in the catchment and burial in a system whereby detrital minerals could either constitute the primary host phase or correlate to ${\mathrm{Hg}}_{\mathrm{T}}$ and so could explain the progressive elevation in ${\mathrm{Hg}}_{\mathrm{T}}$ and ${\mathrm{Hg}}_{\mathrm{AR}}$ observed in Lake Ohrid (Fig. 4).

4.3.2

Last Glaciation (35–12 $\mathrm{ka}$; MIS 3 to MIS 2)

The timing, amplitude, and expression of $\mathrm{Hg}$ signals captured in Lake Prespa and Lake Ohrid change significantly between $\sim$ 35 and 12 $\mathrm{ka}$ (Fig. 7). The largest ${\mathrm{Hg}}_{\mathrm{T}}$ and ${\mathrm{Hg}}_{\mathrm{AR}}$ peaks in Lake Ohrid coincide with the Last Glacial Maximum (LGM) and begin at $\sim$ 35 $\mathrm{ka}$ (Fig. 7). Synchronous enrichments in $\mathrm{K}$, quartz, and $\mathrm{Ti}$ (Francke et al., 2016; Wagner et al., 2019) provide evidence of elevated clastic terrigenous matter inputs and erosion and are consistent with evidence of a significantly less-vegetated catchment (Donders et al., 2021; Sadori et al., 2016). High clastic fluxes into the lake during the LGM could also relate to meltwater runoff from local mountain glaciers (Ribolini et al., 2011), which would transport large volumes of sediment generated by glacial abrasion, quarrying, and plucking (Carrivick and Tweed, 2021; Overeem et al., 2017) into the lake basin. Given that $\mathrm{Hg}$ sequestration in Lake Ohrid appears partially related to the abundance of detrital minerals for much of the record (Figs. 4 and 5), these $\mathrm{Hg}$ peaks could relate to local, climate-driven shifts in landscape structure associated with glaciation during MIS 2 (Figs. 4 and 7).

Alternatively or in addition to these local effects, atmospheric mineral dust concentrations were also up to 20 times higher during the LGM (Simonsen et al., 2019). Mineral dust may be the most important $\mathrm{Hg}$ carrier in ice cores (Jitaru et al., 2009; Vandal et al., 1993), and studies have shown evidence of notable redistribution of terrestrial $\mathrm{Hg}$ during the LGM owing to changes in regional atmospheric dust deposition (de Lacerda et al., 2017; Fadina et al., 2019; Pérez-Rodríguez et al., 2015). However, we see no clear evidence that atmospheric dust played a major (direct) role in the local $\mathrm{Hg}$ cycle in our data. For example, peaks in elemental ratios typically associated with mineral dust deposits (e.g. $\mathrm{Zr}/\mathrm{Ti}$) do not correspond to peaks in ${\mathrm{Hg}}_{\mathrm{T}}$ and/or ${\mathrm{Hg}}_{\mathrm{AR}}$ (Fig. S7 in the Supplement) (Vogel et al., 2010). We also do not see loess-based evidence of elevated aeolian dust fluxes over central Europe and the Balkans during the Last Glacial Maximum (Újvári et al., 2010; Rousseau et al., 2021). Marine sediment records also do not capture measurable changes in Saharan dust influx to the Ionian and Aegean seas, corresponding to pronounced $\mathrm{Hg}$ signals in Lake Ohrid (Fig. S7) (Ehrmann and Schmiedl, 2021). Therefore, we cannot mechanistically link elevated $\mathrm{Hg}$ values during MIS 2 in Lake Ohrid to broad-scale changes in atmospheric dust deposition.

The largest ${\mathrm{Hg}}_{\mathrm{T}}$ and ${\mathrm{Hg}}_{\mathrm{AR}}$ peaks in Lake Prespa occur between 21.3 $\pm$ 1.7 (1$\mathit{\sigma }$ from the Bayesian age model; see Fig. 3) and 17.5 $\pm$ 0.7 $\mathrm{ka}$. These signals do not correspond to a measurable change in host-phase availability (Fig. 5), so it is unlikely that these peaks reflect changes in TOC, TS, and/or $\mathrm{K}$. However, they do coincide with deglaciation of the Pindus and Dinaric mountains (Fig. 7) (Hughes et al., 2023). Geomorphological evidence suggests that glaciers were present across the Prespa–Ohrid region between $\sim$ 26.5 and 15 $\mathrm{ka}$ (Belmecheri et al., 2009; Gromig et al., 2018; Ribolini et al., 2018; Ruszkiczay-Rüdiger et al., 2020) and indeed that periglacial processes created a landscape characterized by intense weathering, erosion, and sediment transport (Hughes and Woodward, 2017; Allard et al., 2021). Glacial meltwaters thus likely constituted a major source of water input to Lake Prespa during the last deglaciation. Glaciers are important sinks for atmospheric $\mathrm{Hg}$ deposited by both dry and wet processes (Durnford and Dastoor, 2011; Zhang et al., 2012), and large quantities of $\mathrm{Hg}$ can accumulate in organic-rich frozen soils (permafrost; Schuster et al., 2018). High proportions of detrital matter within glacial ice, snow, and organic matter facilitate the effective, long-term ($>$ hundreds to thousands of years) retention of atmospheric $\mathrm{Hg}$, meaning that rapid snow/ice melt and permafrost thawing can produce transient “pulses” of $\mathrm{Hg}$ into lakes without a comparable peak in sediment influx (Durnford and Dastoor, 2011; Kohler et al., 2022). This is consistent with the abrupt and short-lived increase in $\mathrm{Hg}$ concentration retained in Lake Prespa between 21.3 and 17.5 ($\pm$ 1.7–0.7, 1$\mathit{\sigma }$) $\mathrm{ka}$, which occurs in the absence of a pronounced change in terrigenous elements (e.g. $\mathrm{Ti}$, $\mathrm{Rb}$) or TS (Figs. 5 and 7).

Lakes Ohrid and Prespa show two other striking differences in $\mathrm{Hg}$ concentration between 35–12 $\mathrm{ka}$ (Fig. 7). First, Lake Prespa does not record a distinct ${\mathrm{Hg}}_{\mathrm{T}}$ or ${\mathrm{Hg}}_{\mathrm{AR}}$ signal during the LGM, and second, Lake Ohrid does not record a distinct ${\mathrm{Hg}}_{\mathrm{T}}$ or ${\mathrm{Hg}}_{\mathrm{AR}}$ signal corresponding to deglaciation. Given their close proximity and environmental similarity, both lakes could be expected to record similar overall signals if the climate-driven processes influencing ${\mathrm{Hg}}_{\mathrm{AR}}$ were broadly similar. One plausible explanation could be a disproportionately large change in Lake Prespa's total volume compared to Lake Ohrid. Increased abundance of small Fragilariaceae and benthic Eolimna submuralis diatom species points to generally low temperatures and lake levels during MIS 2 (Cvetkoska et al., 2015). These conditions are also indicated by elevated concentrations of ice-rafted coarse sand and gravel grains and further suggest persistent ice formation on the lake surface, likely facilitated by the lake's shallow depth (Damaschke et al., 2013; Wagner et al., 2010; Vogel et al., 2010). It is possible that the heightened presence of ice at the peak of glaciation served as a natural barrier between the surface and the sediments of Lake Prespa, effectively slowing the net flux of $\mathrm{Hg}$ into delivery of solutes to the basin. A simultaneous lack of ice cover on Lake Ohrid, linked to greater water depths, could also justify why ${\mathrm{Hg}}_{\mathrm{AR}}$ remained high in this lake during the LGM, as the $\mathrm{Hg}$ influx pathway would be unaffected by ice formation (Fig. 7).

Water volume changes may also have influenced the hydrological connection between lakes Ohrid and Prespa during deglaciation (Cvetkoska et al., 2016; Jovanovska et al., 2016; Leng et al., 2010). Tracer experiments and stable isotope (${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$) analysis suggest that water draining from Lake Prespa accounts for a significant proportion of Lake Ohrid's water inflow alongside precipitation (Matzinger et al., 2006; Wagner et al., 2010; Lacey and Jones, 2018), with high rates of prior calcite precipitation occurring in the connecting karst system (Eftimi et al., 1999; Leng et al., 2010; Matzinger et al., 2006). However, a change to lower ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ of lake water and TIC in both lakes during the last glaciation point to a reduction in the contribution of karst-fed waters to Lake Ohrid (Lacey et al., 2016; Leng et al., 2013). Although it is unlikely that the two hydrological systems became completely decoupled (Belmecheri et al., 2009; Lézine et al., 2010), evidence of permafrost formation at high elevations between 35 and 18 $\mathrm{ka}$ (Oliva et al., 2018) and lower precipitation could be linked to a reduction in karst aquifer activity (Fig. 7). For shallower Lake Prespa, lower precipitation may also have led to a larger reduction in lake volume compared to Lake Ohrid, a decrease in the number (and pressure) of active sinkholes, and subsequently the outflow of water and solutes (e.g. $\mathrm{Hg}$) into Lake Ohrid (Wagner et al., 2014a) – increasing both ${\mathrm{Hg}}_{\mathrm{T}}$ and ${\mathrm{Hg}}_{\mathrm{AR}}$. Together, the collective impact of disproportionately large, climate-driven reductions in water level could explain why rates of $\mathrm{Hg}$ accumulation were significantly higher in Lake Prespa during deglaciation compared to the LGM. Glacial meltwaters would elevate the net $\mathrm{Hg}$ input compared to the LGM, and reduced ice cover would permit a more direct pathway for $\mathrm{Hg}$ to be delivered into the basin, both processes becoming effective, while underground permafrost continued to limit the intra-basin exchange of water and solutes.

Neither Lake Ohrid nor Lake Prespa shows large changes in $\mathrm{Hg}$ concentration or accumulation during the Oldest (17.5–14.5 $\mathrm{ka}$) and Younger (12.9–11.7 $\mathrm{ka}$) Dryas. Both lakes contain clear evidence of an abrupt return to glacial conditions during this time. Lake Prespa sediments record shifts in tree pollen and diatom assemblages, alluding to a net reduction in local winter temperatures and moisture availability (Aufgebauer et al., 2012a; Panagiotopoulos et al., 2013; Cvetkoska et al., 2014), and high uranium (${}^{\mathrm{234}}\mathrm{U}/{}^{\mathrm{238}}\mathrm{U}$) activity ratios, low tree pollen percentages, and low TIC concentrations in Lake Ohrid also pertain to intense hillslope erosion owing to a more open catchment structure (Francke et al., 2019b; Lézine et al., 2010). Geomorphological evidence also pertains to local glacier stabilization (Gromig et al., 2018; Ribolini et al., 2018; Ruszkiczay-Rüdiger et al., 2020) (Fig. 7). Nonetheless, we suggest these events may have been too (a) short-lived and/or (b) climatically mild to produce a similarly distinct response in the terrestrial $\mathrm{Hg}$ cycle as the processes operating during, and immediately following, the LGM, potentially explaining the lack of an associated sedimentary $\mathrm{Hg}$ signal.

4.3.3

Holocene (12–0 $\mathrm{ka}$; MIS 1)

The timing and amplitude of ${\mathrm{Hg}}_{\mathrm{T}}$ and ${\mathrm{Hg}}_{\mathrm{AR}}$ signals recorded in Lake Prespa and Lake Ohrid sediments are noticeably different during the Holocene (MIS 1). Between 12 $\pm$ 0.5 and 3 $\pm$ 0.2 $\mathrm{ka}$, Lake Prespa captures a series of large peaks in ${\mathrm{Hg}}_{\mathrm{T}}$ and ${\mathrm{Hg}}_{\mathrm{AR}}$, corresponding to high TOC and TIC and indicative of elevated productivity, higher rates of organic material preservation, and limited mixing (Fig. 5). Conversely, ${\mathrm{Hg}}_{\mathrm{T}}$ and ${\mathrm{Hg}}_{\mathrm{AR}}$ show a progressive decline in Lake Ohrid during MIS 1, despite coeval increases in TOC and TIC (Fig. 6). These observations suggest that for most of the Holocene, $\mathrm{Hg}$ fluxes into the two lakes were largely decoupled, likely due to differences in catchment and basin dynamics, which impacted the rate of $\mathrm{Hg}$ delivery to (and burial in) the lakes.

Divergent $\mathrm{Hg}$ signals in Lake Ohrid and Lake Prespa during this time may be linked to heightened wildfire frequency and/or intensity. Wildfires have the capacity to (in)directly release $\mathrm{Hg}$ from vegetation and/or through associated changes in soil erosion. Proxy evidence alludes to interglacial conditions characterized by heightened seasonality, marked by very warm, dry summers coupled with wet, mild winters; an overall increase in the prevalence of deciduous tree species (Cvetkoska et al., 2014; Panagiotopoulos, 2013); and also an increase in macro- and microcharcoal concentrations in Lake Prespa (Fig. 7; Panagiotopoulos et al. 2013). Large wildfires would have a broadly regional-scale impact, which, given the close proximity of our two lakes, could theoretically produce a measurable $\mathrm{Hg}$ signal in both systems. However, more frequent and/or intense regional fires could also yield measurably different sedimentary $\mathrm{Hg}$ signals by their capacity to (1) enhance surface runoff without a corresponding increase in erosion and effectively reduce transport of catchment-sourced, mineral-hosted $\mathrm{Hg}$ (Mataix-Solera et al., 2011; Shakesby, 2011); (2) enhance downstream transport of $\mathrm{Hg}$ released from burned soils and bound to fine and coarse particulate matter (Burke et al., 2010; Takenaka et al., 2021); and/or (3) release large quantities of $\mathrm{Hg}$ into the atmosphere following biomass combustion (Howard et al., 2019; Melendez-Perez et al., 2014; Roshan and Biswas, 2023). All three combine to generate impacts that may vary in significance owing to lake-specific differences in sedimentation, accumulation, and flux of materials to/from the lake.

An increase in wildfire activity also corresponds to a period of intensifying human influence in the region, predominantly in the form of land use change, agriculture, and animal husbandry (Cvetkoska et al., 2014; Masi et al., 2018; Panagiotopoulos et al., 2013; Rothacker et al., 2018; Thienemann et al., 2017; Wagner et al., 2009). Widespread mineral resource exploitation and metalworking on the Balkan Peninsula have been recorded as early as $\sim$ 8 $\mathrm{ka}$ (Gajić-Kvašèev et al., 2012; Longman et al., 2018; Radivojević and Roberts, 2021; Schotsmans et al., 2022), and release of detrital $\mathrm{Hg}$ during cinnabar ore extraction and use of $\mathrm{Hg}$ in gold extraction (amalgamation) have been linked to pronounced $\mathrm{Hg}$ contamination in modern sedimentary units in the region (Covelli et al., 2001; Fitzgerald and Lamborg, 2013). Directly quantifying the influence of (hydro)climate- versus human-driven impacts on sedimentary $\mathrm{Hg}$ records presents a major challenge as these factors are interdependent. Nonetheless, these factors could produce a more measurable effect in lake systems with heightened sensitivity to changes in water, nutrients, and pollutant fluxes. This could explain why large $\mathrm{Hg}$ signals are observed in Lake Prespa between $\sim$ 12 and 3 $\mathrm{ka}$ but not in Lake Ohrid: Lake Prespa is shallow relative to its surface area (Fig. 2), meaning that relatively small oscillations in pollutant influxes can lead to appreciable changes in lake geochemistry (Cvetkoska et al., 2015; Matzinger et al., 2006).

Decoupling of the two $\mathrm{Hg}$ records effectively disappeared $\sim$ 3 $\mathrm{kyr}$ ago, during which both lakes showed a sharp and pronounced rise in ${\mathrm{Hg}}_{\mathrm{T}}$ and ${\mathrm{Hg}}_{\mathrm{AR}}$ (Fig. 7). Several lines of evidence point to human activity as the primary cause. On a local scale, a rapid increase in the biological productivity (eutrophication) of Lake Prespa since $\sim$ 1.6 ($\pm$ 0.06) ka alludes to greater disturbance of catchment soils by agricultural practices and eventually the use of inorganic compounds such as pesticides and fertilizers (Aufgebauer et al., 2012; Cvetkoska et al., 2014; Krstić et al., 2012; Leng et al., 2013). Signals observed in Fig. 7 may thus be a product of human-induced changes in organic or minerogenic material flux: each facilitating more efficient delivery of catchment-sourced $\mathrm{Hg}$ (Fitzgerald et al., 2005) and possibly also stimulating microbial $\mathrm{Hg}$ methylation within the sediment (Soerensen et al., 2016). On a broader scale, peaks in ${\mathrm{Hg}}_{\mathrm{T}}$ and ${\mathrm{Hg}}_{\mathrm{AR}}$ correspond to a sustained rise in European and/or global $\mathrm{Hg}$ emissions, owing to increased deforestation, fossil fuel extraction and combustion, and intentional use of $\mathrm{Hg}$ for resource extraction/production (Outridge et al., 2018; United Nations Environment Programme, 2018). An increasing number of sedimentary archives record $\mathrm{Hg}$ enrichments as early as $\sim$ 3 $\mathrm{kyr}$ ago (Biskaborn et al., 2021; Guédron et al., 2019; Li et al., 2020; Pan et al., 2020). The emergence of simultaneous ${\mathrm{Hg}}_{\mathrm{T}}$ and ${\mathrm{Hg}}_{\mathrm{AR}}$ peaks in lakes Ohrid and Prespa following $\sim$ 3 $\mathrm{ka}$ underscores the magnitude and global distribution of this change in $\mathrm{Hg}$ sources and emissions (Fig. 7) and points to a rise in $\mathrm{Hg}$ fluxes between 3 and 0 $\mathrm{ka}$ that was distinct enough to effectively overwhelm previously dominant natural drivers of $\mathrm{Hg}$ variability.

4.4 Key differences and implications

The magnitude and expression of $\mathrm{Hg}$ signals recorded in Lake Prespa and Lake Ohrid are different in three aspects. First, the extent to which different host phases can (or cannot) explain time-varying patterns in $\mathrm{Hg}$ concentration differs between the two lakes. Although only a limited fraction of $\mathrm{Hg}$ variability in either record can be explained by availability of any single host phase, the low degree of covariance that we do observe points to organic material playing the most significant role as a $\mathrm{Hg}$ host in Lake Prespa. In contrast, $\mathrm{Hg}$ correlates most strongly with detrital minerals in Lake Ohrid over the same period (0–90 $\mathrm{ka}$) (Fig. 4). The second difference is visible during the last glaciation ($\sim$ 35–12 $\mathrm{ka}$): in Lake Ohrid $\mathrm{Hg}$ concentrations peak during the LGM (35.8–12 $\mathrm{ka}$), whereas Lake Prespa captures transient, high-amplitude peaks during deglaciation, starting $\sim$ 15 $\mathrm{kyr}$ later (Fig. 7). The third difference is visible during the Holocene. The largest signals in the entire Lake Prespa record are observed between $\sim$ 8 and 0 $\mathrm{ka}$, whereas $\mathrm{Hg}$ concentrations do not increase in Lake Ohrid until $\sim$ 2 $\mathrm{ka}$. These observations raise the following question (Fig. 2): for two lakes located in such close geographical proximity and having experienced similar climate conditions, what may have caused such pronounced differences from $\sim$ 35 ka?

Differences in bathymetric structure may offer a plausible explanation. For example, the largest changes in the amplitude and frequency of peaks in ${\mathrm{Hg}}_{\mathrm{T}}$ and ${\mathrm{Hg}}_{\mathrm{AR}}$ are exhibited by Lake Prespa (Fig. 7): a shallow basin that contains $>$ 90 % less water than Lake Ohrid, despite only a $\sim$ 30 % difference in surface area (Wagner et al., 2010). Increased distance from the lake margin to the core site in Lake Ohrid would mean distribution of material over a greater total area and thus more time for net $\mathrm{Hg}$ loss to occur by evasion from the water surface (Cooke et al., 2020), removal of water (and suspended material) via riverine outlets (Bishop et al., 2020), or processes taking place within the water column (Frieling et al., 2023) prior to burial. Therefore, preservation of a measurable $\mathrm{Hg}$ signal in a deep lake (e.g. Lake Ohrid) would require notably larger influx of $\mathrm{Hg}$, and this sedimentary signal would also likely be significantly smaller than the equivalent “dose” delivered to a smaller and/or shallower lake (e.g. Lake Prespa). Coupled with evidence of high-amplitude fluctuations in lake water ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}$ ($\pm$ 6 ‰) (Leng et al., 2010) and lake level (Cvetkoska et al., 2015, 2016) corresponding to pronounced $\mathrm{Hg}$ variability in Lake Prespa, but not in Lake Ohrid (Fig. 7), our data suggest that smaller, shallower lakes may be particularly sensitive recorders of transient changes in $\mathrm{Hg}$ fluxes.

Divergent bathymetric structures are also linked to distinct differences in biological composition and nutrient availability in lakes Ohrid and Prespa. The deep ($\sim$ 240 $\mathrm{m}$) waters of Lake Ohrid host a highly oligotrophic (nutrient-poor) environment characterized by low levels of biological productivity and a high abundance of planktonic diatom species (e.g. Cyclotella) (Cvetkoska et al., 2021). Conversely, Lake Prespa's shallower ($\sim$ 14 $\mathrm{m}$) waters host a dominantly mesotrophic (nutrient-rich) system in which benthic and planktonic diatom species are present in equal abundance (Jovanovska et al., 2016; Cvetkoska et al., 2016) and allude to moderate/high biological productivity (Leng et al., 2013). Productivity is a potentially important factor influencing the $\mathrm{Hg}$ composition of lake sediment: high productivity typically favours higher concentrations of algal biomass, allowing for more effective $\mathrm{Hg}$ scavenging by organic particles and export to the sediment (Biester et al., 2018; Soerensen et al., 2016; Hermanns et al., 2013). While the overall signal will remain dominated by $\mathrm{Hg}$ availability, broad-scale differences in productivity between lakes Prespa and Ohrid through time could provide an additional explanation for the disparate expression of recorded $\mathrm{Hg}$ signals (Sect. 4.1), with notably higher productivity in the shallower Lake Prespa further increasing its sensitivity to changes in nutrient status, erosion, and hydrology.

Local differences in $\mathrm{Hg}$ emission by neotectonic activity may also have contributed to the divergent $\mathrm{Hg}$ signals, owing to differences in the host rock geology, tectonic instability, and mechanical stress regimes of faults surrounding the two basins (Hoffmann et al., 2010; Lindhorst et al., 2015). However, the significance of these differences cannot be fully assessed in the absence of direct $\mathrm{Hg}$ emission measurements (see Text S8 in the Supplement).

The two records presented here highlight that $\mathrm{Hg}$ cycling in lacustrine environments is distinct from open marine systems. In marine systems, $\mathrm{Hg}$ fluxes can be broadly modulated by large-scale continental sediment (Fadina et al., 2019; Figueiredo et al., 2022; Kita et al., 2016) and/or atmospheric inputs (Chede et al., 2022), and $\mathrm{Hg}$ burial flux ultimately becomes more closely related to host-phase availability. Conversely, both Lake Prespa and Lake Ohrid highlight how the local basin and catchment characteristics both exert key control over the delivery of $\mathrm{Hg}$ to lacustrine sediments and suggest that differences in $\mathrm{Hg}$ cycling between geographically proximal basins could occur as a function of diverse physical, hydrological, and biological properties.

Our observations highlight that multi-millennial lacustrine $\mathrm{Hg}$ records allow a different perspective of the $\mathrm{Hg}$ cycle compared to marine records and, for example, may be used to infer how local, regional, and global climatic conditions could have altered processes important to the terrestrial $\mathrm{Hg}$ cycle. Because lacustrine records are much better suited to recording smaller-scale processes, it is also clear that extrapolating the (non-marine) $\mathrm{Hg}$ cycle response from a single lacustrine $\mathrm{Hg}$ record is challenging. For example, a single-core approach could produce a large degree of uncertainty owing to variable sediment focussing and catchment-sourced influx of organic and inorganic materials (Blais and Kalff, 1995; Engstrom and Rose, 2013; Engstrom and Wright, 1984). A valuable next step would be to apply a source-to-sink approach within a well-known lacustrine catchment: to assess the extent to which $\mathrm{Hg}$ sedimentation is spatially heterogeneous within a lacustrine system and whether multiple cores extracted from different locations within the same basin would yield markedly different $\mathrm{Hg}$ trends. Intra-basin heterogeneity in $\mathrm{Hg}$ sources, reactions, and transformations could also be examined through measurement of stable $\mathrm{Hg}$ isotopes, particularly in millennium-scale sedimentary records where the nature of these processes may change through time (Blum et al., 2014; Jiskra et al., 2022; Kurz et al., 2019). Work of this nature would make great strides toward assessing how representative of variability in the local $\mathrm{Hg}$ cycle a single, in this case, lake core is, and whether intra-basin fluctuations in sedimentation, resuspension, and erosion could translate to measurable changes in sedimentary $\mathrm{Hg}$ burial.

Past changes in environmental $\mathrm{Hg}$ availability inferred from sedimentary records have typically been examined (and presented) by normalizing $\mathrm{Hg}$ to a dominant host phase, often taken as organic matter (Fadina et al., 2019; Figueiredo et al., 2020; Grasby et al., 2019; Kita et al., 2016; Percival et al., 2015). However, availability of organic matter or other host phases that scavenge $\mathrm{Hg}$ here appear to represent just one of several processes governing $\mathrm{Hg}$ burial in lacustrine systems, and this process is very likely systematically less significant compared to marine records in lieu of changes in catchment and basin processes such as erosion, nutrient status, and hydrology (Outridge et al., 2019). Outside pre-industrial times (or periods without an overwhelming global $\mathrm{Hg}$ cycle perturbation, such as during large igneous province (LIP) formation; Grasby et al., 2019), a single common process/mechanism is therefore unlikely to produce a unanimous stratigraphic signal across all lakes or even for two adjacent lakes as shown in this study.

5 Conclusions

To better understand the local and regional impacts of climate, vegetation, and catchment characteristics on lacustrine $\mathrm{Hg}$ records, we present two new high-resolution $\mathrm{Hg}$ records for the last $\sim$ 90 $\mathrm{kyr}$ from Lake Prespa and Lake Ohrid. The two records show some similarities but also distinct differences in the strength of the relationships among $\mathrm{Hg}$, TOC, TS, and detrital minerals ($\mathrm{K}$), with only a relatively small proportion of $\mathrm{Hg}$ variability attributable to host-phase availability in each record. Our findings provide three valuable insights. First, the local sedimentary environment does influence $\mathrm{Hg}$ burial. Covariance with host phases accounts for a limited proportion of the observed variability, suggesting that many of the ${\mathrm{Hg}}_{\mathrm{T}}$ and ${\mathrm{Hg}}_{\mathrm{AR}}$ signals recorded in Lake Prespa and Lake Ohrid reflect net $\mathrm{Hg}$ input to the two lakes across timescales ranging from decades to multiple millennia. Second, $\mathrm{Hg}$ signals can reflect changes in (and also differences between) catchment hydrology and structure. Despite their proximity, the magnitude and expression of the recorded signals are considerably different between Lake Prespa and Lake Ohrid, suggesting these inputs changed relative to the sedimentary setting and in response to changing interactions between the two systems. Finally, regional-scale climate variability can measurably affect the $\mathrm{Hg}$ signals retained in lake sediments: both Lake Prespa and Lake Ohrid show changes in $\mathrm{Hg}$ concentration and accumulation, corresponding to glacial (Late Pleistocene) and interglacial (Holocene) climate conditions. It follows that local, regional, or global changes in climate or hydrological cycling capable of affecting mineral soils, (peri-)glacial features, or fire regime in the lake catchment could all impact $\mathrm{Hg}$ fluxes. These findings prompt further examination of how orbital-scale climate variability ($>$ 10${}^{\mathrm{3}}$-year timescales) may influence the terrestrial $\mathrm{Hg}$ cycle, not only to better resolve processes acting on single lacustrine and terrestrial successions, but also to identify which of these (local) processes could hold relevance for $\mathrm{Hg}$ cycling on a global scale.

Code availability

Details of materials and software used in generation of the Co1215 chronology are provided in the text and in the Supplement. This includes the rBacon software, which is available online (https://rdrr.io/cran/rbacon/).

Data availability

Data generated as part of this study will be available in the Supplement. Pre-existing data for lakes Prespa and Ohrid used in this study are available either as supplementary datasets in the cited literature or via the PANGAEA online repository (https://www.pangaea.de/, last access: 15 September 2023). Material from cores Co1215 and 5045-1 are stored partly at the University of Cologne and at the MARUM core repository in Bremen, Germany, and are available on request.

The supplement related to this article is available online at: https://doi.org/10.5194/bg-21-531-2024-supplement.

Author contributions

ARP, IMF, JF, TAM, DMP, and SAR conceptualized the study. TAM acquired the funding. ARP performed Hg analyses for Co1215 and 5045-1 with assistance from IMF. JL and BW acquired and provided core material for Hg analysis. IMF generated the new age model for Co1215. JL, AF, BW, TT, and KP provided additional geochemical data for both cores analysed in this study. ARP wrote the original manuscript draft, and all authors contributed to review and editing.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Special issue statement

This article is part of the special issue “Land surface–atmosphere interactions – from the microbial to the global scale”. It is not associated with a conference.

Acknowledgements

Alice R. Paine, Isabel M. Fendley, Joost Frieling, and Tamsin A. Mather acknowledge funding from the European Research Council Consolidator Grant V-ECHO (ERC-2018-COG-8187 17-V-ECHO). Alice R. Paine thanks David Thomas and Mona Edwards (School of Geography, Oxford) for logistical assistance with sample transfer and storage. Konstantinos Panagiotopoulos acknowledges funding from the German Research Foundation (DFG; grant PA 2664/4-1). All authors thank members of the Scientific Collaboration on Past Speciation Conditions in Lake Ohrid (SCOPSCO) and the CRC 806 Our Way to Europe – Culture-Environment Interaction and Human Mobility in the Late Quaternary projects for their efforts in producing the Lake Ohrid and Lake Prespa sediment successions and making the data available for scientific use.

Financial support

This research has been supported by a European Research Council Consolidator Grant (V-ECHO, grant no. ERC-2018-COG-8187 17-V-ECHO) and funding from the German Research Foundation (DFG; grant no. 2664/4-1).

Review statement

This paper was edited by Semeena Valiyaveetil Shamsudheen and reviewed by three anonymous referees.