A global compilation of diatom silica oxygen isotope records from lake sediment – trends and implications for climate reconstruction

Oxygen isotopes in biogenic silica (δ¹⁸O_BSi) from lake sediments allow for quantitative reconstruction of past hydroclimate and proxy-model comparison in terrestrial environments. The signals of individual records have been attributed to different factors, such as air temperature (T_air), atmospheric circulation patterns, hydrological changes, and lake evaporation. While every lake has its own local set of drivers of δ¹⁸O variability, here we explore the extent to which regional or even global signals emerge from a series of paleoenvironmental records. This study provides a comprehensive compilation and combined statistical evaluation of the existing lake sediment δ¹⁸O_BSi records, largely missing in other summary publications (i.e. PAGES network). For this purpose, we have identified and compiled 71 down-core records published to date and complemented these datasets with additional lake basin parameters (e.g. lake water residence time and catchment size) to best characterize the signal properties. Records feature widely different temporal coverage and resolution, ranging from decadal-scale records covering the past 150 years to records with multi-millennial-scale resolution spanning glacial–interglacial cycles. The best coverage in number of records (N = 37) and data points (N = 2112) is available for Northern Hemispheric (NH) extratropical regions throughout the Holocene (roughly corresponding to Marine Isotope Stage 1; MIS 1). To address the different variabilities and temporal offsets, records were brought to a common temporal resolution by binning and subsequently filtered for hydrologically open lakes with lake water residence times < 100 years. For mid- to high-latitude (> 45° N) lakes, we find common δ¹⁸O_BSi patterns among the lake records during both the Holocene and Common Era (CE). These include maxima and minima corresponding to known climate episodes, such as the Holocene Thermal Maximum (HTM), Neoglacial Cooling, Medieval Climate Anomaly (MCA) and the Little Ice Age (LIA). These patterns are in line with long-term air temperature changes supported by previously published climate reconstructions from other archives, as well as Holocene summer insolation changes. In conclusion, oxygen isotope records from NH extratropical lake sediments feature a common climate signal at centennial (for CE) and millennial (for Holocene) timescales despite stemming from different lakes in different geographic locations and hence constitute a valuable proxy for past climate reconstructions.

1 Introduction

1.1 Scientific background

Oxygen isotopes are ubiquitous within the global water cycle and are among the best (hydro)climate proxies worldwide as a result of their potential quantitative interpretability. The most abundant isotope ¹⁶O and the rarer isotope ¹⁸O are subject to fractionation during water-phase transformation and transport processes (Fig. 1). As a result, the relative abundance of these isotopes varies across space, time, and reservoirs. The IAEA's standard Vienna Standard Mean Ocean Water (VSMOW) serves as a baseline, closely resembling the isotopic composition of modern seawater. Relative abundance of ¹⁸O with regard to ¹⁶O is expressed relative to the VSMOW standard as δ¹⁸O and given in units of per mill.

Figure 1 Schematic overview of select reservoirs and processes within the global water cycle influencing the δ¹⁸O_BSi signal of records. (a) Overview of the global water cycle with reservoirs and processes. (b) Depleting and enriching factors influencing δ¹⁸O values in precipitation, lake water, and biogenic silica. Corresponding fractionation processes are indicated on the right-hand side.

Measuring δ¹⁸O in water, vapour, snow, and ice allows for quantifying hydroclimatic processes in the global water cycle and has been used to this end for many decades (Dansgaard, 1964; Konecky et al., 2020). In addition to its use in present-day hydrology, δ¹⁸O in environmental archives also serves as a powerful proxy for past hydrology and, in turn, climate. As such, δ¹⁸O is a crucial and quantitative tool for both reconstructions of past climate and for comparison of proxy data to climate model outputs.

Oxygen isotope records are available from different paleoenvironmental archives, e.g. in glaciers and ice sheets (Andersen et al., 2004), permafrost (Meyer et al., 2015b), and marine (Spielhagen and Mackensen, 2021) and lacustrine (Leng, 2006) environments, as well as different materials, such as carbonates (Kwiecien et al., 2014), silicates (Leng, 2006), biomarkers (Lasher et al., 2017), cellulose (Wolfe et al., 2007), glacier ice (Andersen et al., 2004), and ground ice (Opel et al., 2017). Therefore, oxygen isotopes offer the possibility for directly comparing data from different paleoenvironmental archives.

The predictability of isotope fractionation processes in the water cycle due to clear physical constraints allows for using δ¹⁸O not only as a quantitative tool for past (hydro)climate reconstructions but also for implementing δ¹⁸O in global climate models (e.g. Danek et al., 2021). However, challenges persist in understanding and comparing data and model outputs, partly due to complex signal formation in environmental archives (Danek et al., 2021).

Lake sediments are prominent and widespread archives in terrestrial environments (Leng, 2006; Biskaborn et al., 2016; Subetto et al., 2017). Like glaciers and ice caps, lakes may record the isotopic signature of past precipitation (Shemesh et al., 2001a). The signal formation in lakes, however, is fundamentally different from that in glaciers and ice caps, with the lake waterbody acting as a buffer of the enclosed (hydro)climate signal, influenced by factors such as catchment hydrology, lake ontogeny, and sediment accumulation rates (Bittner et al., 2021). While biogenic silica is also found in phytoliths, sponge spicules, and chrysophytes, this study focuses exclusively on diatoms. Diatoms are microscopic algae that synthesize their frustules from lake water and the silica dissolved therein, thereby preserving the isotope signal of past lake water. The oxygen isotope composition of the biogenic silica of diatom frustules (δ¹⁸O_BSi) buried in lake sediment therefore provides valuable insight on past lake water and, in turn, precipitation.

Naturally, the interplay of hydrological and sedimentological processes limits the time span and resolution over which environmental information can be obtained from a lake system, and this varies from lake to lake. Most commonly, δ¹⁸O_BSi has been applied to lake sediments devoid of carbonates. These lakes are especially common in high-altitude and high-latitude regions (Leng and Barker, 2006).

1.2 Controls on δ¹⁸O_BSi

Due to the complex signal formation in δ¹⁸O_BSi records, most authors refer to a combination of factors for interpreting any given individual record. A schematic overview of some of the main processes is provided in Fig. 1 and by Leng and Barker (2006). Lake water temperature (T_lake) imparts a direct effect on δ¹⁸O_BSi due to temperature-dependent fractionation during biosynthesis of biogenic silica. Temperature-dependent fractionation coefficients were initially determined for marine environments (Leclerc and Labeyrie, 1987; Matheney and Knauth, 1989) and later estimated for lacustrine sediments as ca. −0.2 ‰ per °C during biosynthesis (Dodd and Sharp, 2010; Moschen et al., 2005). However, the applicability of these calibrations to sedimentary records is debated due to possible disequilibrium between diatom frustules and lake water, as well as taphonomic and diagenetic processes (Leng and Barker, 2006; Ryves et al., 2020; Smith et al., 2016; Tyler et al., 2017).

More commonly, changes in δ¹⁸O_BSi are attributed to δ¹⁸O changes in lake water (δ¹⁸O_lake) and thereby, with the lake as a buffer, to the isotopic composition of precipitation (δ¹⁸O_{precipitation}). Locally, δ¹⁸O_{precipitation} is primarily influenced by air temperature (T_air) with a fractionation of 0.7 ‰ per °C on global average (Dansgaard, 1964). On a broader scale, parameters such as precipitation moisture source and air mass trajectory patterns are important, and changes in δ¹⁸O_BSi records may be indicative of shifting atmospheric circulation patterns (Bailey et al., 2015) that affect the local ratio of precipitation to evaporation ($P/E$). The parameter $P/E$ refers to the balance between the amount of precipitation and evaporation within a lake's catchment. Typically, water loss through evaporation leads to higher δ¹⁸O_lake values due to evaporative enrichment and therefore higher δ¹⁸O_BSi. Precipitation amount and seasonality (e.g. rain vs. snow) will also strongly impact the isotope signal of the lake water (Meyer et al., 2022). Disentangling the effects of precipitation and evaporation in δ¹⁸O_BSi records is thus a major challenge, and therefore some authors commonly refer to $P/E$ changes in conjunction with other factors (Hernandez et al., 2013; Rosqvist et al., 2013; Broadman et al., 2020a).

Shemesh et al. (2001a) have attributed the δ¹⁸O_BSi signal changes of Lake 850 in Swedish Lapland to different contributions of main source regions of precipitation, namely the Atlantic and Arctic oceans. Indeed, for this location, changes in atmospheric circulation patterns imply changes in both T_air and humidity.

Additionally, δ¹⁸O_{precipitation} may also be linked to shifts in the seasonality of precipitation (Swann et al., 2010; Kostrova et al., 2016; Harding et al., 2020), which is linked to T_air change, but also atmospheric circulation patterns (Shemesh et al., 2001a; Jones et al., 2004; Leng et al., 2005; Rosqvist et al., 2013) and, on longer timescales, insolation (Swann et al., 2010; Kostrova et al., 2019, 2021). T_air as often inferred from δ¹⁸O_BSi can, therefore, be a parameter with local, regional, and global significance.

Insolation is not an environmental parameter that can be inferred from δ¹⁸O_BSi records. It is, however, commonly used as an explanatory variable and underlying driver of Holocene climate change. Insolation can be calculated for any given place throughout time (Laskar et al. 2004), which allows for direct comparison with δ¹⁸O_BSi time series. Insolation primarily influences T_air and, indirectly, the $P/E$ balance and atmospheric circulation patterns.

Moreover, hydrological processes within a lake's catchment may substantially impact δ¹⁸O_lake and hence δ¹⁸O_BSi. Some of these processes are closely linked to climate and precipitation patterns described above. Most notable are variations in the amount of snowmelt (Mackay et al., 2013; Rosqvist et al., 2013; Meyer et al., 2022) and glacial influx (Meyer et al., 2015a) to a lake. Glaciers or snow fields in the hinterland may provide more melt (with lower or depleted δ¹⁸O) in warm phases and may counter the influence of warming on δ¹⁸O_{precipitation} (Meyer et al., 2015a).

Other hydrological processes, however, are only indirectly linked to climate. These include the formation and closure of outflows to the lake which change the lake's hydrological setting (Hernandez et al., 2008; Vyse et al., 2020; Bittner et al., 2021). Associated changes in the mass balance of a lake, and thus the $P/E$-balance, may lead to a different δ¹⁸O_BSi signal. Within the lake itself, long-term diagenetic effects on diatoms include recrystallization and incorporation of heavier ¹⁸O from silanol bonds into the Si–O–Si crystal lattice – a process that may increase the original δ¹⁸O_BSi signal (cf. Akse et al., 2022; Fig. 1). However, despite notable recrystallization, Chapligin et al. (2012b) found no significant diagenetic effect in records spanning the past 250 kyr.

To interpret a given δ¹⁸O_BSi record, information about the present hydrology is thus important as it provides insight on which processes influence δ¹⁸O_lake and, in turn, δ¹⁸O_BSi. Naturally, as the hydrology of a given lake can only offer snapshots of present-day constraints, caution has to be applied when extrapolating these into the past. The isotopic composition of three components is of interest for assessing the signal properties of δ¹⁸O_BSi records: (1) δ¹⁸O_lake, (2) δ¹⁸O of inflows to the lake, and (3) δ¹⁸O_{precipitation}. The lake water composition can provide insight on whether lake water is isotopically homogenous, both spatially and at different depths of the lake basin. Offsets of δ¹⁸O_lake from the Local Meteoric Water Line (LMWL) give hints as to whether or not lake evaporation may have had a significant effect on lake water isotopic composition. δ¹⁸O_lake can also be used in conjunction with T_lake and the most recent δ¹⁸O_BSi value for testing the temperature-dependent fractionation and, in turn, the applicability of the proxy for paleoclimate reconstructions (cf. Meyer et al. 2015a). Complementary isotope measurements of lake inflows (i.e. rivers, streams) provide further information on specific hydrological constraints. For larger lakes and catchments, the different inflows and their effect on the lake's isotopic mass balance have been used for inferring changing hydrology or precipitation regimes in different parts of the catchment (e.g. Mackay et al. 2011). Additionally, while δ¹⁸O_{precipitation} data are critical for constraining modern climate–isotope relationships, they are often not available and only few monitoring studies of δ¹⁸O_{precipitation} and δ¹⁸O_lake in combination with δ¹⁸O_BSi exist (Kostrova et al., 2019; Hernandez et al., 2010). A possible alternative can be derived from the GNIP database (IAEA/WMO, 2022) and spatial interpolations thereof, e.g. the Online Isotopes in Precipitation Calculator (Bowen, 2022).

1.3 Sample preparation and measurement

Extracting pure diatom valves (frustules) from lake sediments is a complex process and depends on a sufficiently high diatom concentration in the sediment. Various cleaning procedures for diatom extraction exist, which are mainly based on the use of H₂O₂ for removing organic matter from the sediment and HCl for the removal of carbonate (Chapligin et al., 2012a). Detrital components are often separated through centrifugation in heavy liquid solution (i.e. sodium polytungstate) at variable specific densities (Morley et al., 2004; Chapligin et al., 2012a). Different procedures in sample preparation, i.e. different temperatures, may result in offsets of measured δ¹⁸O_BSi (Swann et al., 2010; Chapligin et al., 2012b; Tyler et al., 2017). Impurities due to incomplete removal of detrital components may also affect the measured δ¹⁸O_BSi (Lamb et al., 2005; Chapligin et al., 2012b). To account for this, the purity of samples is usually assessed prior to measurement by visual inspection or direct measurement of the amount of contaminants of a given sample (Bailey et al., 2018; Broadman et al., 2020a). Fewer studies determine the isotopic composition of detrital contaminants to apply a correction to the measured δ¹⁸O_BSi values (Brewer et al., 2008; Mackay et al., 2011; Wilson et al., 2014a; Bittner et al., 2021; Kostrova et al., 2021). Consequently, caution has to be applied when comparing absolute values of different individual δ¹⁸O_BSi records, as offsets between individual records may result from both different preparation and measurement techniques, as well as regional differences in δ¹⁸O_lake. The potential and challenges associated with data stemming from different preparation and measurement techniques have already been addressed (Chapligin et al., 2011; Mackay et al., 2011).

1.4 Aim of this work

A comprehensive compilation and assessment of the lacustrine δ¹⁸O_BSi records published to date is missing. Due to the crucial role of lake and catchment hydrology in signal formation, such a compilation needs to include individual lake basin parameters, such as lake volume and water residence time (t_res), to provide reliable constraints for interpreting and comparing the individual records. Proxy data compilations have already addressed a lack of standardized metadata, data availability, and data uniformity of paleo-data (Pfalz et al., 2021). However, such compilations generally do not include δ¹⁸O_BSi records and their associated lake basin parameters (Kaufman et al., 2020) or instead they have a limited temporal focus (Konecky et al., 2020). Consequently, no study has yet empirically linked the signal properties of δ¹⁸O_BSi records and lake basin parameters in a harmonized dataset.

To overcome these gaps, this paper aims at providing a comprehensive compilation and combined statistical evaluation of all lake sediment δ¹⁸O_BSi records published to date (between 1998 and 2022). We accomplish these objectives by means of the following working steps: (1) collecting available lake sediment δ¹⁸O_BSi records published to date, (2) complementing these records with the individual lake basin parameters, (3) assessing the signal properties of δ¹⁸O_BSi records with regard to lake basin parameters, and (4) identifying spatio-temporal patterns and trends in the δ¹⁸O_BSi signals. This effort leads to a better understanding of the general constraints for interpreting lake sediment δ¹⁸O_BSi records, whilst making these records more readily usable for future studies including proxy-model comparisons. It is hence a contribution to bridge the gap between modelling and isotope geochemistry approaches in paleoclimate science.

2 Methods

This study follows a three-stage approach: the first stage (data acquisition) comprises identifying lacustrine δ¹⁸O_BSi datasets and publications published to date. A second stage includes acquiring the actual δ¹⁸O_BSi datasets and archiving them in a standardized format. Where appropriate, the identified published datasets were arranged into longer continuous records. This includes records from the same sediment core published in different papers and records from different sediment cores from the same lake but similar locations. The third and final stage (record analysis) interprets the lacustrine δ¹⁸O_BSi isotope records with respect to their hydrological and geographical constraints in order to identify possible common global or hemispheric trends and signal properties of all isotope records, as well as smaller spatially or temporally constrained datasets.

Identification and acquisition of the datasets and publications in this work followed an additive approach, i.e. thorough literature survey. We have chosen this approach due to the limited number of laboratories and working groups worldwide measuring δ¹⁸O_BSi. In this study, we focused exclusively on down-core records from lacustrine sediments. Identified publications and datasets were entered into a uniform metadata table giving one entry to each publication and each dataset. If more than one publication was written about the same sediment core, each of these publications were given a separate entry. Likewise, if one publication presents data from more than one sediment core, each of the sediment cores were given a separate entry to the database.

For each of these entries, metadata on coring procedure, hydrological setting, and chronology were supplemented. Data were extracted from the corresponding publication(s), from public repositories, or directly from the authors. Hydrological parameters such as catchment area, average depth, and t_res were additionally obtained from the HydroLAKES database (Messager et al., 2016) and stored as separate variables. For a detailed description on how the values in the HydroLAKES database were created, we refer to Messager et al. (2016). This procedure was necessary because original publications do not always specify all of these parameters and because parameters supplied in original publications may not be consistent. For further analysis and for linking the individual isotope records with lake basin parameters, the HydroLAKES database entries were given preference. Original publications were only used for lakes where no data were available from the HydroLAKES database. The parameter “maximum water depth” was always taken from the original publications since this parameter is not provided by the HydroLAKES database. While the geo-statistical approach of the HydroLAKES database may not provide the most precise values for individual lake basin characteristics, it ensures comparability among the different lake basins analysed in this study. Therefore, and because some parameters (i.e. catchment size) are not essential to the discussion in this paper, we have chosen this approach favouring consistency over more data points. It also prevents potential issues arising from different authors providing conflicting values for the same lake basin (i.e. Lake Baikal, Lake El'gygytgyn).

Isotope datasets for each of the entries were archived in separate tables that specify depth, age, and measured isotope values. Datasets were taken directly from the respective publication, where possible. If datasets were not available, within the original publications or their supplements, public repositories (e.g. http://www.pangaea.de, last access: August 2022; http://ncei.noaa.gov, last access: August 2022) were searched. In instances where datasets were unavailable from repositories, the lead authors were contacted directly. In case of published data unavailable from repositories or authors, plots of the original publications were digitized if possible (Rietti-Shati et al., 1998; Hu and Shemesh, 2003). Where digitizing plots was not feasible, records were excluded from further consideration; this was the case for Chondrogianni et al. (2004) and Hu et al. (2003).

Chronologies were adapted from the original publications, where available, and stored in cal yr BP format (relative to 1950 CE). Sample ages given in different formats (e.g. b2K or CE) were converted to cal yr BP. This procedure also applies to radiocarbon-based chronologies. Chronologies were not recalculated because the effect of different ¹⁴C calibrations and different age model approaches is presumably minor with regard to the aim of this study. This is especially true when considering the fact that high-resolution datasets (annual to decadal resolution) are the exception among all datasets identified. The error introduced by the different age models is thus considered minor when compared to the resolution of the datasets. For applications requiring more precise chronologies and better comparability of datasets, we provide the radiocarbon measurements of respective datasets as well. This enables future users to create tailor-made chronologies if needed. Datasets without chronologies were stored using depth notation only.

For further analysis, datasets were partly regrouped and combined to generate longer continuous δ¹⁸O_BSi datasets, herein referred to as “records”. This comprised treating datasets with data from several coring sites or outcrops from the same lake as one single record, in agreement with the author's original interpretation (Quesada et al., 2015; Swann et al., 2018). Data stemming from the same sediment core but published in different papers were combined into single continuous records. This applies to records from Lake Kotokel (Kostrova et al., 2013a, b, 2014, 2016) and Lago Chungará (Hernandez et al., 2008, 2010, 2011, 2013). Likewise, data from different cores presented in different publications but stemming from identical or reasonably similar coring sites were combined to single records. This applies to data from Lake Nar (Dean et al., 2018) and Vuolep Allakasjaure (Rosqvist et al., 2004; Jonsson et al., 2010). Lakes consisting of several basins were generally treated as one, such as Lake Baikal and several smaller lakes.

For trend analysis and inter-comparison, we focused primarily on the Holocene and Common Era climate using records from Northern Hemisphere (NH) extratropical lakes (45–90° N) due to there being the highest spatial–temporal coverage of data with a significant number of records, which makes an interpretation departing from case studies feasible.

In a first step, records were binned to 1 kyr intervals (for the Holocene) and 200 year intervals (for the CE), respectively. These intervals were chosen based on the temporal resolutions of the original records in order to ensure continuous binned records with no empty bins. Choosing higher resolutions would have resulted in many records having empty bins with no data points, thus leading to different temporal resolutions even after binning. The binning was done by calculating the mean value of all samples within the respective age interval for each individual record. Data points with ages < 150 yr BP were excluded from the Holocene analyses to eliminate any possible effects of recent warming, though these data were retained for analysing the Common Era. To assess the timing of Holocene maxima and minima, records covering fewer than 10 bins were discarded. To evaluate Holocene and Common Era trends in δ¹⁸O_BSi, records covering less than seven bins were excluded from the analysis (i.e. < 7.0 and 1.4 kyr, respectively). Obviously, the missing bins in both Holocene and Common Era compilations might have an effect on the mean value of the individual records but not on the overall trend.

In order to eliminate offsets between individual binned records, the remaining records were standardized by subtracting their respective means. After this mean removal, two subsequent filtering steps were performed in order to exclude records most prone to secondary lake evaporation. First, a subset consisting only of records from open lakes was created, discarding records from semi-closed lakes, closed lakes, and paleo-lakes. In a second step, this open-lakes subset was filtered by discarding all records with t_res> 100 years. This threshold was chosen in accordance with previous works classifying lakes and their isotopic signals with respect to hydrologic setting and t_res (Leng and Marshall, 2004).

After applying these filters, combined trends for geographical regions (NH, Eurasia, North America) were evaluated by calculating the mean of all records for each bin, in each respective region.

Figure 2 Temporal coverage by records and data points during MIS 1–MIS 5. Boundaries between MIS according to Lisiecki and Raymo (2005): MIS 1–2: 14 kyr; MIS 2–3: 29 kyr; MIS 3–4: 57 kyr; MIS 4–5: 71 kyr; MIS 5–6: 130 kyr. (a) Total number of data points available during each Marine Isotope Stage. (b) Temporal coverage of records from young to old.

3 Results

3.1 Published datasets

Following a thorough literature survey, a total number of 71 published down-core datasets of δ¹⁸O_BSi from 64 sites have been identified. Since the first published lacustrine δ¹⁸O_BSi studies in 1998 (Rietti-Shati et al., 1998; Shemesh and Peteet, 1998), there has been a growing research interest and number of publications (Fig. A1). Only recently has there been a stagnation of the number of records published. A detailed overview of the identified publications and the records presented therein is provided in Table A1.

Most publications present the δ¹⁸O_BSi data as time series, as is usual in paleoenvironmental studies. The chronologies rely on a wide variety of dating methods (Fig. A2), with ¹⁴C (used on 46 records) being by far the most frequent. For sub-recent time periods and shorter timescales, ²¹⁰Pb and ¹³⁷Cs are also used extensively (on 21 and 12 records, respectively). Globally or regionally correlated time markers such as tephra layers (Heyng et al., 2015) or high-resolution methods such as varve counts (Rozanski et al., 2010) are used relatively sparsely as they are not available for all time periods and/or lake basins.

3.2 Combined records

We combined separate records from the same lake sites to longer, continuous composite records, giving a total of 53 combined δ¹⁸O_BSi records. Six of these are from paleo-lakes (three with or without chronology), while 47 records stem from extant lakes (42 with and five without chronology). An overview of these combined records is provided in Table A2, complemented by metadata of the records and the corresponding lake basins. Each individual record received a number (#X; Table A2, column 1; Table A1, last column) which is used consistently throughout the text.

Figure 3 World map depicting the locations of δ¹⁸O_BSi records with chronology. Note that one location may correspond to more than one record. © Esri 2022.

The temporal coverage of these combined records (Fig. 2) varies from centennial-scale at Laguna Zacapu (#50, Leng et al., 2005) and a specific study at Lake Baikal (#51, Swann et al., 2018), to glacial–interglacial cycles, e.g. 250 kyr at Lake El'gygytgyn (#3, Chapligin et al., 2012b). Most records (N = 32) cover the last several thousand years up to about 10 kyr BP, while only a few records cover several tens of thousands of years (e.g. #8 to #10). Likewise, the time periods covered do vary either due to the availability of lake sediments, diatom abundances in a lake and/or scientific foci of the research groups. By far the highest number of records is available for the Holocene (0–11.7 kyr BP; N = 37) and especially for the last 2 kyr BP (N = 48). Further back in time, Marine Isotope Stages (MISs) were used for assessing data coverage throughout time. Boundaries between MISs are according to Lisiecki and Raymo (2005). MIS 2 is still covered by 19 records, whereas towards MIS 3 (N = 5) and beyond fewer records have been generated. However, single records cover even older time periods in MIS 11 (at Lake Baikal; #2) and beyond (MIS G1/104, #1), outlining the applicability of this proxy across a wide range of timescales. A paleolake from the Baringo–Bogoria basin (#1) is the oldest lacustrine δ¹⁸O_BSi record and has been dated to the onset of MIS 104 in the late Pliocene (Wilson et al., 2014b).

Regarding the number of published records by continent (Fig. 3), there is a clear concentration on Asia and Europe with 21 and 15 records published, respectively. In Asia, there is a strong regional focus on Siberia, whereas most other parts of the continent are yet to be investigated using lake sediment δ¹⁸O_BSi. A regional focus also occurs in Africa and South America. While a sizeable number of records stemming from these continents have been published (N = 10 and N = 11, respectively), there are pronounced regional foci in the East African Rift and the Andes, respectively. Another focus region is Alaska, where most (N = 6) of the published North American records (N = 7) are located. δ¹⁸O_BSi work has also been carried out in lakes from New Zealand (#42) and South Georgia (#19).

In summary, the distribution of available δ¹⁸O_BSi data displays a bias towards the mid- to high latitudes of the Northern Hemisphere. Additionally, this distribution is also indicative of the site accessibility and (regional) research foci of the working groups. This has a rather pronounced effect on the spatial distribution of available δ¹⁸O_BSi records. The geographical distribution is also indicative of cold regions devoid of carbonates where biogenic silica is the most promising archive to obtain oxygen isotope records from lake sediments. Moreover, high northern latitudes have by far the highest abundance of total area of waterbodies (Downing and Duarte, 2009). Yet, it also accounts for the relatively lower number of northern records extending beyond MIS 1, as vast ice sheets covered much of the NH high latitudes (e.g. Patton et al., 2016).

Figure 4 (a) Depiction of corresponding lake and catchment sizes of the records compiled within this study. One lake may correspond to more than one record. Note that records from paleo-lakes and lakes with incomplete information on lake and catchment size were not considered for this figure (N = 33). (b) Depiction of corresponding sampling intervals and t_res. Note that records from lakes without information on t_res were not considered for this figure (N = 37).

This pattern also affects the latitudinal distribution of the records (see also Table A2) with records ranging from 54.17° S (#19) to 69° N (#39). Records stemming from low-latitude lakes (particularly in Africa and South America) are often located at high altitudes above 3000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (above sea level; e.g. #11, 17, 18, 21). Lake Simba Tarn (record #34) features the maximum altitude for an individual δ¹⁸O_BSi record of 4959 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$. Altitudes below 1000 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$, however, are most common (N = 36), especially for high-latitude lakes. In total, 12 lakes are located below 100 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (e.g. #29, 35). Many of these low-altitude lakes are located in maritime locations in immediate proximity to the coast (e.g. #19, 30) or have even had marine intrusion stages in the past (a postglacial transgression at Nettilling Lake; #41). Extremely continental environments are tackled in central and eastern Siberia and the Lake Baikal region (e.g. #4, 8, 9, 10, 25).

3.3 Spatial coverage and resolution of combined records

Varying spatial resolution among the compiled isotope records is linked to the fact that lakes and their corresponding catchments are not points in space but areas. Consequently, lake water and the related δ¹⁸O_BSi signal represent a spatial average integrated over the catchment size of each respective lake. Lake surface areas range from small ponds of < 1 km² (e.g. #19, 11, 23, 33, 34) to some of the largest lakes in the world, including Lake Malawi and Lake Baikal covering 29 544 and 31 968 km² , respectively (Fig. 4a). Large and voluminous lakes are the exception, however, and most of the lakes represented here are < 10 km² in size (N = 25). Note that such parameters are not applicable to paleolakes. While some authors provide estimates on paleolake extent, these figures are not consistent with the values determined for present lakes and are therefore not included in the evaluation of the dataset. Catchment sizes vary by several orders of magnitude as well. For small lakes, catchments are often < 10 km² (N = 13; e.g. #22, 33, 37), whereas the largest lakes, Lake Malawi and Lake Baikal, feature catchments of 128 727 and 569 176 km², respectively (Fig. 4a). Therefore, these lakes integrate the environmental signal over a large and potentially diverse hydrological region. While most of the lakes compiled in this study are primarily fed by surface runoff and/or precipitation (according to the original publications of records), groundwater influx may also play a pivotal role by introducing a large memory effect of past precipitation due to the generally long residence times of aquifers. This may have an impact on records especially when looking at short timescales. Groundwater input, however, is usually not accounted for and is thus beyond the scope of our study. In summary, both lake and catchment sizes vary by several orders of magnitude (Fig. 4a) among the sites with existing δ¹⁸O_BSi records and span from local signals to regional averages. However, most of the records (N = 18) stem from lakes with catchments < 100 km², suggesting rather local signals. While single local signals represent small areas, different local signals may well correlate on continental and hemispheric scales.

3.4 Temporal resolution of combined records

Temporal resolution of the records and the resulting signal properties are determined by both the lake basin itself (i.e. accumulation rates and preservation of diatom silica) and the sampling routine applied to the sediment core (i.e. the sampled intervals as well as the thickness of individual samples). Data on the thickness of single sediment samples taken from the cores are, however, scarce in the original publications. The same applies for information on the time interval or number of years represented by a single sediment sample. Consequently this parameter and its potential effect on the records' signal properties could not be investigated in detail in this paper. We therefore focus on t_res and sampling frequency of a core as a means of characterizing the records and their temporal resolution.

The t_res is closely linked to the size of lake basins and varies accordingly. It ranges from several weeks for small lakes (e.g. #49, 50) to centuries (219 and 375 years for Lake Malawi and Lake Baikal, respectively). On the whole, t_res of sub-annual to annual scale are most common (N = 17) among the lakes considered in this study (Fig. 4b). It should be noted that longer t_res leads to averaging of the signal over a longer time period. Likewise, different t_res may also have an effect on absolute δ¹⁸O_BSi values and variabilities. Lakes with very long t_res (t_res> 100 years) are more susceptible to lake evaporation, and may effectively behave like closed-system lakes with regard to the δ¹⁸O_BSi signal even if they are hydrologically open (i.e. have outflows). These are still referred to as hydrologically open in this work and hydrological settings and t_res are addressed separately in the discussion.

The sampling frequency also varies by several orders of magnitude (Fig. 4b). It stretches from annual to multi-millennial timescales. Consequently, the signal properties and the recorded climatic and/or hydrological forcings of the identified records can be expected to vary accordingly. Most records, however, plot in the top left half of graph of Fig. 4b, which indicates that the temporal offset between two sediment samples exceeds t_res. This suggests that the sampling routine is the limiting factor of the temporal resolution of these records. Only three records from Lake Baikal plot in the lower-right half of Fig. 4b, suggesting t_res to be the dominant factor in determining the record's resolution in these cases (#4, #20, #51), at least with respect to inflow-related changes.

Figure 5 Common Era Northern Hemispheric records (45–90° N) compiled in this study. (a) Original records, (b) data binned to 200-year intervals only showing records covering at least seven bins, (c) binned records with mean of individual records removed, (d) data filtered for records from open lakes only, (e) data filtered for records from lakes with t_res < 100 yr, and (f) the NH trend calculated as mean of all records in each bin. Shadings show 1 and 2 standard deviations, respectively. (g) North American and (h) Eurasian trend.

We thus conclude that the sampling resolution is the main factor determining the temporal resolution of most records, with t_res acting as an additional smoothing mechanism. A lake with a centennial-scale t_res can, therefore, display centennial-scale changes in climate and hydrology to their full amplitude. Decadal-scale changes, on the other hand, can be expected to be attenuated by an order of magnitude in a lake with centennial-scale t_res. When t_res is lower than the sampling frequency, the hydroclimatic amplitude is also not fully captured (Richter and Turekian, 1993). It should be noted that temporal resolution is non-uniform across the records, with generally higher resolution for more recent time intervals (Figs. 2 and 5a) that – in addition to a sampling bias – may in part reflect increasing lake sediment compaction with depth and time. Based on t_res and sampling resolutions of the records, comparing records on a centennial or millennial scale is the most promising approach for assessing common patterns.

4 Discussion

4.1 Common Era climate

For analysis of Common Era climate, only data points from 0 CE until today (corresponding to 1950 to −73 yr BP) were selected. The Common Era subset of records includes only records stemming from sites north of 45° N (Fig. 5a, N = 19) and comprises 460 data points. The records display considerable offsets with δ¹⁸O_BSi values ranging from +19 ‰ (#34) to +33 ‰ (#35) VSMOW. These offsets can be linked to their individual environmental setting, e.g. latitudinal and continentality effects, as well as to their potential to be prone for lake evaporation. Amplitudes in δ¹⁸O_BSi of CE records vary from 0.4 ‰ (#24) to ca. 9 ‰ VSMOW (#35). This variability might correspond to their different hydrological settings. Closed lakes (such as Sunken Island Lake; #26) have a tendency towards higher δ¹⁸O_BSi amplitude due to evaporative enrichment of the lake water (Broadman et al., 2020b).

Binned records (Fig. 5b, N = 14) do not display any common pattern and feature δ¹⁸O_BSi values ranging from +20 ‰ to +31 ‰ VSMOW. However, binning of data points and exclusion of shorter records reduces the δ¹⁸O_BSi amplitudes of individual records to less than 5 ‰. Offsets between individual records may be linked to site-specific characteristics of individual lakes and catchments, as well as latitudinal and continentality effects.

After mean removal (i.e. standardization; Fig. 5c), the CE data suggest a common trend with higher δ¹⁸O_BSi values at 1900 yr BP and a general decrease until the present, with the lowest values occurring in the last 400 years (two bins). Some records also show maxima between 500 and 1100 yr BP and minima at ca. 1300 yr BP. The subsets for hydrologically open lakes (N = 13, Fig. 5d) and hydrologically open lakes with t_res < 100 yr (N = 9, Fig. 5e) show a similar pattern.

The combined CE records (Fig. 5f) comprise eight records and show a general δ¹⁸O_BSi decrease over the last 2 kyr amounting to ca. 2 ‰ VSMOW (slope: ∼ 1 ‰ kyr⁻¹). The δ¹⁸O_BSi amplitude within individual bins, however, might exceed this number, most notably at 100, 1300, and 1900 yr BP, which show the highest amplitudes of up to 5 ‰ for the CE. This is in line with differences between the timing of δ¹⁸O_BSi maxima and minima between North American and Eurasian records. At centennial scale, these differences might also be linked to dating uncertainties. Moreover, the record of Lake Bolshoye Shchuchye (#16) features a much larger amplitude. Its exceptional δ¹⁸O_BSi variability has been attributed to variations in snow and snowmelt in the lake's catchment and therefore does represents a different kind of precipitation-based signal compared to most other records. The combined trend features its highest δ¹⁸O_BSi values at 1900 yr BP, followed by a decrease leading to a relative minimum at 1100 and 1300 yr BP. This phase is followed by a second δ¹⁸O_BSi peak at 700 and 900 yr BP, after which a decrease can be observed. The decrease at 300 and 500 yr BP features the least deviation between individual records. The most recent bin again shows a much larger δ¹⁸O_BSi variability. While most records display the lowest values at this time, Lake Kotokel (#8) for instance shows an increase compared to the previous bins. This increase may or may not be related to recent warming, and due to the complex hydrology of the lake even lake evaporation cannot be ruled out. Most of the records, however, do not indicate a recent δ¹⁸O_BSi maximum that might be indicative of recent warming.

It should be noted that there is bias of data points towards the older end of the most recent bin. An absence of recent warming in the data therefore might be linked to this bias, and hence a reaction of lake systems to recent warming cannot be ruled out. As Eurasian records constitute the majority of records (N = 7), it is similar to the NH trend (Fig. 5h).

The trends calculated for the last 2 kyr (2K) for the NH (Eurasian) δ¹⁸O_BSi stacks are high with −0.64 ‰ kyr⁻¹ (−0.85 ‰ kyr⁻¹) for 10 bins, with even higher slopes of −1.6 ‰ kyr⁻¹ (−1.7 ‰ kyr⁻¹) for bins 6–10 (1 to 2 kyr) and −1.4 ‰ kyr⁻¹ (−1.8 ‰ kyr⁻¹) for bins 1–5 (corresponding to the last millennium). The negative slope is interrupted by three consecutive bins (4–6; 0.6–1.2 kyr BP) plus bin 8 (1.4–1.6 kyr BP) with higher δ¹⁸O_BSi than the previous bin. For the last 2 kyr BP, the other five bins (1–3, 7, and 9) show a negative sign for the NH and Eurasian δ¹⁸O_BSi 2K reconstructions.

North American records (N = 2) do not show a consistently decreasing trend (Fig. 5g) but show slightly higher values at 1.7 and 1.9 kyr BP compared to the most recent bins (0.1 and 0.3 kyr BP, respectively). They do, however, show lower δ¹⁸O_BSi values between 0.9 and 1.3 kyr BP, followed by a δ¹⁸O_BSi maximum at 500 yr BP. As there are only two records available after filtering, caution has to be applied in interpreting this pattern.

For the North American δ¹⁸O_BSi 2K reconstruction, a lower overall gradient is observed (−0.12 ‰ kyr⁻¹), which shows also a steep slope for bins 6–10 with −1.5 ‰ kyr⁻¹. In contrast to the NH and Eurasian δ¹⁸O_BSi stack, bins 1–5 show a slight decrease of −0.4 ‰ kyr⁻¹ only. This leads to slightly shifted minima and maxima between δ¹⁸O_BSi reconstructions for the last 2 millennia: NH and Eurasian reconstructions have their absolute maxima in bin 10 or at 1.8–2.0 kyr (absolute minima: bin 1; 0–0.2 kyr) with intermediate minima at 1.2–1.4 kyr BP (bin 7) and maxima between 0.6–1.0 kyr BP (bins 4 and 5). In contrast, for the North American reconstruction, the absolute δ¹⁸O_BSi minimum (maximum) is at 1.2–1.4 kyr BP or bin 7 (0.4–0.8 kyr BP (bins 3 and 4), whereas the early maximum (bin 10) and late minimum (bin 1) are less pronounced. In summary, we observe an overall decreasing trend in δ¹⁸O_BSi for the last 2 millennia for the NH, Eurasian, and North American stacks, which is accelerated in the first millennium.

Similar patterns among Eurasian and possibly NH δ¹⁸O_BSi records do suggest a common NH signal throughout the CE. The observed maxima and minima correspond to previously described climatic events, notably the Roman Warm Period from 2.0 kyr BP to about 1.6 kyr BP (Ljungqvist, 2010), the Dark Ages Cold Period from 1.6 to 1.2 kyr BP (Büntgen et al., 2016; Helama et al., 2017), the Medieval Climatic Anomaly (MCA) from 1.2 to 0.8 kyr BP (Bradley et al., 2003; Mann et al., 2009), and the Little Ice Age from 0.7 to 0.1 kyr BP (Matthews and Briffa, 2005; Wanner et al., 2022). In Fig. 5f, the two main cold phases (Dark Ages Cold Period and Little Ice Age) are clearly represented. The good accordance of our data with these previously described warm and cold phases suggests that δ¹⁸O_BSi records presented are influenced by T_air change, either directly or via other parameters linked to T_air. Recent research has rejected the global nature of these climatic events and suggested that they are regionally constrained (Neukom et al., 2019). However, the accordance of our data – when assessed at centennial scale – with these climatic events is expected because they were initially described for North America and Europe.

Figure 6 Holocene Northern Hemispheric records (45–90° N) compiled in this study. (a) Original records, (b) data binned to 1 kyr intervals only showing records covering at least 7 bins, (c) binned records with mean of individual records removed, (d) data filtered for records from open lakes only, (e) data filtered for records from lakes with t_res < 100 yr, and (f) the NH trend calculated as mean of all records in each bin. Shadings show 1 and 2 standard deviations, respectively. (g) North American and (h) Eurasian trend.

4.2 Holocene (corresponding to MIS 1) climate

To facilitate comparison with other Holocene reconstructions, for deriving Holocene trends only 12 kyr (in 1 kyr bins) are considered, and the boundary between MIS 1 and MIS 2 is set to 12 kyr. Thus, the Holocene NH subset of records (covering the past ca. 12 kyr) consists of 21 records that display a considerable range in mean δ¹⁸O_BSi values from +20 ‰ (#3) to +35 ‰ (#31) VSMOW (Fig. 6a) that reflect local site-specific difference. For example, Lake El'gygytgyn (#3) is located in a continental, high-latitude environment with little secondary evaporation (Chapligin et al., 2012b), whereas Lake Ladoga (#31) is situated in a less continental, lower-latitude setting with a complex hydrological history of changing inflows and outflows (Kostrova et al., 2019).

Amplitudes of the individual Holocene records vary from less than 5 ‰ (#3) to ca. 10 ‰ VSMOW (#20, #35). Lake El'gygytgyn (#3) is a voluminous, deep hydrologically open lake with little lake evaporation, whereas Heart Lake (#35) is a much smaller, less voluminous lake with a short t_res (2 weeks). Lake Baikal (#20) is a deep voluminous lake like Lake El'gygytgyn, but it has a centennial-scale t_res making it effectively closed and thus potentially subject to secondary evaporation. Consequently, differences in their isotopic variability might correspond to their different hydrological settings. Closed lakes (such as Sunken Island Lake; #26) have a tendency towards higher δ¹⁸O_BSi values that are likely due lake evaporation leading to a more enriched isotope signature of lake water (Broadman et al., 2020b; Hernandez et al. 2008).

Binned records (Fig. 6b, N = 14) do not display any common patterns and feature δ¹⁸O_BSi values ranging from +20 ‰ to +35 ‰ VSMOW, similar to the original records. However, binning of data points and exclusion of shorter records results in smoothed amplitudes of less than 5 ‰ VSMOW. The binned records exhibit Early Holocene δ¹⁸O_BSi values (up to 2.5 ‰) higher than the Holocene mean, whereas late Holocene values are generally up to 2 ‰ below the Holocene mean. However, some records do not follow this general pattern after mean removal (Fig. 6c). This could be indicative of the different hydrological settings and t_res which have a substantial impact on the recorded signal.

Filtering for hydrologically open lakes (Fig. 6d, N = 13) and t_res < 100 yr (Fig. 6e, N = 10) displays a clear decreasing δ¹⁸O_BSi trend of 2.5 ‰ throughout the Holocene. However, differences between individual records do exist, and most notably the absolute δ¹⁸O_BSi maxima of individual records do not occur at the same time. While some records feature maxima at the beginning of the Holocene, other records peak between 5 and 8 kyr BP. This pattern is in line with Holocene T_air reconstructions which have found spatial differences of the timing of the Holocene thermal maximum (Kaufman et al., 2004; Renssen et al., 2009).

Based on this subset of hydrologically open lakes with t_res < 100 yr, a combined millennial-scale NH trend was calculated. This combined Holocene trend (Fig. 6f) shows a decrease throughout the Holocene amounting to ca. 2 ‰ VSMOW. The absolute maximum is observed for the Early Holocene bin at 11–12 kyr BP. A decrease at the beginning of the Holocene between 12 and 10 kyr BP is followed by a relatively stable middle Holocene until 6 kyr BP and subsequent stronger decrease towards the absolute δ¹⁸O_BSi minimum in the youngest bin (0–1 kyr BP).

When considered separately, North American (N = 2) and Eurasian (N = 8) records show different patterns and have been described and interpreted differently in the case studies published to date. While the former are primarily linked to atmospheric circulation changes in Alaska (Bailey et al., 2018; Broadman et al., 2020b, a), the latter are at least partly interpreted as indicative of T_air and insolation changes (Swann et al., 2010; Mackay et al., 2011; Chapligin et al., 2012b; Kostrova et al., 2013a, 2019, 2021). Eurasian δ¹⁸O_BSi records (N = 8) display a slight decrease from 12 kyr BP to 10 kyr BP and suggest a second maximum at 7 kyr BP followed by a second relatively stable phase. An accelerated decrease starting at ca. 4 kyr BP constitutes a more abrupt change than in the NH trend. The Holocene trend calculated for the NH (Eurasian) δ¹⁸O_BSi stack is −0.19 ‰ kyr⁻¹ (−0.21 ‰ kyr⁻¹) for all 12 bins, with a lower slope of −0.11 ‰ kyr⁻¹ (−0.10 ‰ kyr⁻¹) between 12–7 kyr BP and a much higher slope of −0.39 ‰ kyr⁻¹ (−0.36 ‰ kyr⁻¹) for bins 1–6 (corresponding to 6–0 kyr BP).

Throughout the Holocene, 10 out of 12 bins are lower than the previous one (i.e. negative sign) for the NH δ¹⁸O_BSi reconstruction, except between 9–8 and 7–6 kyr BP (bins 9 and 7), which show an increase in δ¹⁸O_BSi. For the Eurasian δ¹⁸O_BSi stack all bins are lower in δ¹⁸O_BSi than the preceding one (except for bin 7, or 7–6 kyr BP, which shows an increase of +0.2 ‰). The most negative slopes exceeding −0.55‰ kyr⁻¹ are observed between 2–1 kyr BP for both the NH and the Eurasian compilation. For the North American δ¹⁸O_BSi reconstruction, a slightly lower trend is observed for the Holocene (−0.13 ‰ kyr⁻¹, based on 10 bins), which also shows a steeper slope between 5–0 kyr BP (bins 5–1) with −0.22 ‰ kyr⁻¹. In contrast to the NH and Eurasian δ¹⁸O_BSi stack, there is a slight increase of +0.10 ‰ kyr⁻¹ between 10–6 kyr BP (bins 10–6), likely driven by bins 6 and 9, which are the only ones showing a positive sign, whereas all 7 other bins show a decrease compared to the preceding one.

Hence, North American records do not show a consistently decreasing trend throughout the Holocene (Fig. 6g) but instead slightly higher Δ¹⁸O_BSi values in the first half of the Holocene as compared to the second half. Since there are only two North American records fulfilling our selection criteria, it is difficult to meaningfully go beyond the existing case studies, outlining the necessity for further research in this region. In summary, we observe a mostly continuous decrease in δ¹⁸O_BSi throughout the Holocene for the NH, Eurasian and North American stacks, which is accelerated in the second half of the Holocene. This phenomenon has been previously described as Neoglacial Cooling (e.g. McKay et al., 2018).

The regional differences between North America and Eurasia are also manifested in the timing of absolute δ¹⁸O_BSi minima and maxima of individual records (Fig. 7). The spatial pattern of the Holocene maxima of binned records shows a different timing of maxima for different regions (Fig. 7). Eastern Eurasian sites feature a pronounced early Holocene maximum (at 12 kyr BP), whereas sites to the west of Eurasia show a tendency towards rather a middle Holocene maximum, around 8–6 kyr BP. This suggests that the double maxima of the Eurasian trend (Fig. 6h) is at least partly caused by the regional differences over Eurasia. There are also individual records, however, which do show two peaks within the Holocene, e.g. Lake Bolshoye Shchuchye from the polar Ural Mountains (#16). Records from Alaska feature later δ¹⁸O_BSi maxima (5 and 8 kyr BP, respectively).

Figure 7 Timing of δ¹⁸O_BSi maxima and minima of binned records in the Holocene. Filter criteria are the same as indicated in Fig. 10, and only records covering at least 10 bins have been considered. © Esri 2022.

The timing of Holocene minima is rather consistent across northern Eurasia, where all records reaching their minimum in the last 1 kyr BP or between 2 and 1 kyr BP (Fig. 7). While the records generally follow the same decreasing long-term δ¹⁸O_BSi trend, they feature either early or late Holocene minima. Since data points of 1850 CE and younger have been removed from the dataset prior to analysis to exclude the industrial era, a potential effect of recent climate change is not covered in this subset of data. As with the Holocene maximum, sites in Alaska also differ with regard to the Holocene minimum. They show a tendency towards absolute minima earlier than Eurasian records (between 3 and 4 kyr BP, respectively). Despite the scarcity of the records fulfilling the selection criteria, this suggests a different behaviour of the regions (North America and Eurasia), though we also acknowledge that one of the two North American sites (#35, Heart Lake) is located in the central northern Pacific Ocean and has a different climatic history.

The African (N = 7) and South American records (N = 3) do not fulfil the criteria set for the NH and Eurasian stacks due to either hydrological or temporal constraints (e.g. too few data points). However, all three South American records (see Table A2; #18, 28, 32) show an early Holocene maximum around 10 kyr BP. Two out of three South American δ¹⁸O_BSi records show a decrease over the Holocene, whereas the third record displays no clear trend.

Due to marked differences in hydrology among sites, an African δ¹⁸O_BSi stack has not been calculated either. Instead, individual records (N = 5) have been compiled and binned separately for the African continent (Fig. A5a and b). We observe the largest δ¹⁸O_BSi variability within one individual record (δ¹⁸O_BSi range of ca. 20 ‰, #21) as well as highest absolute δ¹⁸O_BSi (with values up to +45 ‰, #33). This is related to the very different settings both in altitude (Table A2) and hydrological characteristics (Fig. A5d). Overall, the first half of the Holocene is characterized by slightly lower Δ¹⁸O_BSi values than the second half. As this observation is less obvious for open lakes (Fig. A5d), a bias linked to widely different hydrological settings of the respective lakes (see Table A2) has to be assumed, which renders it difficult to disentangle the drivers of δ¹⁸O_BSi. This further underscores the importance of hydrologically constrained records to infer a common (climate) signal (see Figs. 5a and 6a).

Figure 8 Combined trend of δ¹⁸O_BSi compared to other climate reconstructions. (a) Greenland ice sheet temperature reconstruction. (b) Multi-proxy temperature reconstruction (30–90° N). (c) Global stack of δ¹⁸O from benthic foraminifera. (d) Temperature reconstruction (45–90° N) based on pollen data. (e) Temperature anomalies of mean surface temperatures (multi-method ensemble). (f) Insolation anomaly for June and December. (g) Combined NH trend from δ¹⁸O_BSi records (this work). (h) Combined Eurasian and North American trends from δ¹⁸O_BSi records (this work).

4.3 Combined Holocene trend in the hemispheric context

The combined trends of δ¹⁸O_BSi for the NH, Eurasia, and North America are shown in comparison to other NH proxy records in Fig. 8. As individual δ¹⁸O_BSi records are commonly discussed with respect to insolation, we also include June and December insolation curves for 60° N (Laskar et al., 2004, Fig. 8f). As shown in Fig. 8g and h, Eurasian and NH combined trends show a similarity with the June insolation as all records feature a decreasing trend throughout the Holocene. However, insolation decreases steadily after an early Holocene maximum, whereas combined NH and Eurasian δ¹⁸O_BSi trends feature a stable early Holocene and a second peak at 7 kyr BP. This presumably relates to the geographical distribution of the timing of δ¹⁸O_BSi maxima and minima as discussed above. Therefore, eastern Eurasian records are most in line with June insolation, which has been regarded as a proxy of summer T_air (#3, 24). Records stemming from sites further west are less directly correlated with June insolation. A striking difference is the accelerated decrease of δ¹⁸O_BSi after 4 kyr BP in the NH and Eurasian time series, co-incident with Neoglacial Cooling, which is not mirrored in the insolation curve. The relatively stable insolation during this time suggests that the accelerated decrease visible in δ¹⁸O_BSi records must be driven by other factors.

Conversely, December insolation (Fig. 8f) shows an anticorrelation with the δ¹⁸O_BSi records. However, as both the absolute values and the changes in December insolation are an order of magnitude lower than those of June insolation, a decisive influence of December insolation on δ¹⁸O_BSi records can be ruled out. This is in good agreement with previous works, claiming that the δ¹⁸O_BSi proxy yields a summer-dominated signal (Shemesh et al., 2001a; Kostrova et al., 2021), supported by the fact that most biogenic production is in late spring–early summer (particularly relevant for NH high-latitude short-residence-time systems). It has to be noted that the records displayed still stem from different latitudes, and the insolation patterns are not identical at all these sites. A decrease in summer insolation, however, is generally the case for high-latitude regions throughout the Holocene.

Temperature reconstructions by Vinther et al. (2009) using δ¹⁸O data from ice cores at Agassiz ice cap and Greenland (between 65 and 80° N) show a pronounced T_air increase of about 5 °C at the beginning of the Holocene until 10 kyr BP, followed by a stable phase until 7 kyr BP (Fig. 8a). While the NH and Eurasian δ¹⁸O_BSi trends do not feature this increase during the early Holocene, they do show a relatively stable phase from 10 to 7 kyr BP. After 7 kyr BP, both the NH δ¹⁸O_BSi trend and the NH T_air reconstruction feature a decrease of ca. 2 °C and 1.5 ‰ VSMOW, respectively. As the present-day global average T_air-dependent fractionation in precipitation amounts to 0.695 ‰ (°C)⁻¹ (Dansgaard, 1964), a 1.5 ‰ δ¹⁸O decrease would thus correspond to a 2 °C cooling, the same T_air change as found by Vinther et al. (2009). This good agreement of δ¹⁸O_BSi combined trends (NH and Eurasia) and the T_air reconstructions by Vinther et al. (2009) suggests summer T_air to have a major impact on δ¹⁸O_BSi records on millennial timescales.

Multi-proxy-based temperature reconstructions by Marcott et al. (2013) (regional stack 30 to 90° N, Fig. 8b) also show a temperature increase until 10 kyr BP, followed by a stable phase until ca. 7 kyr BP and a decrease of ca. 2 °C thereafter. This temperature reconstruction includes a much broader geographical focus and different land-based and marine proxies (e.g. pollen, chironomids, ice). Thus, the stable phase from 10 to 7 kyr BP and a maximum around 7 kyr BP are shared features of the δ¹⁸O_BSi trends of this work, the ice core reconstruction from Vinther et al. (2009) and the temperature reconstructions by Marcott et al. (2013). However, the early Holocene maximum of the δ¹⁸O_BSi trends is not represented in reconstructions by Vinther et al. (2009) and Marcott et al. (2013). This discrepancy might reflect different regional biases of Marcott et al. (2013) and this work, which is also in line with the different timing of Holocene maxima in the δ¹⁸O_BSi data (Fig. 8). Given the discrepancy between temperature reconstructions and δ¹⁸O_BSi records in the early Holocene, it is likely that the influence of factors other than T_air – such as glacial retreat and ice melt during deglaciation – is especially pronounced in the early Holocene. The LR04 stack of the δ¹⁸O of benthic foraminifera (Lisiecki and Raymo, 2005) features continuously decreasing δ¹⁸O throughout the Holocene, which is in accordance with the combined Eurasian δ¹⁸O_BSi trend (Fig. 8c).

Pollen records may help to further investigate this issue as they stem from terrestrial environments but are also archived in lake sediments. Moreover, they offer a similar temporal and geographical focus compared to our records, extending to high latitudes. For this comparison we have used the Holocene pollen data compilation by Herzschuh et al. (2021, 2022, 2023a, b) extracting data from 864 records from sites north of 45° N. Changes in annual T_air and July T_air relative to modern values were retrieved from the pollen dataset. Annual and July T_air show similar patterns with a pronounced early Holocene increase, followed by middle Holocene maximum and a less pronounced decrease until the present day (Fig. 8d). Both the middle Holocene maximum and the subsequent decrease in T_air are in good agreement with the combined Eurasian δ¹⁸O_BSi trend (Fig. 8d and h). This further supports a substantial influence of T_air on the δ¹⁸O_BSi signal. However, the early Holocene maximum of δ¹⁸O_BSi in our combined NH and Eurasian trends is not reflected by pollen-based reconstructions. Caution has to be applied when using pollen-based climate reconstructions for comparison because vegetation changes may lag millennia behind changes in climatic variables (Herzschuh et al., 2016). In summary, the similarity of our combined δ¹⁸O_BSi trends with NH temperature reconstructions and insolation data allows for deducing a clear link to summer T_air for the Holocene records. Therefore, our δ¹⁸O_BSi record contributes to the understanding of the Holocene climate history.

A significant discrepancy between Holocene cooling deduced from proxy reconstructions (e.g. Marcott et al., 2013; Kaufman et al., 2020; Fig. 8b and e) and Holocene warming simulated in climate models has been called the Holocene “temperature conundrum” (Liu et al., 2014), which was revisited in Wanner (2021), Kaufman et al. (2020) and Kaufman and Broadman (2023). This discrepancy has been attributed to uncertainties in both proxy reconstructions and climate models. One major aspect is the seasonal bias in organic-based proxy records (such as pollen and diatoms) towards summer (Liu et al, 2014). As the Holocene displays opposite trends in summer and winter insolation at 60° N (Fig. 8e, Laskar et al., 2004), high latitudes provide an optimal setting for testing the seasonality aspect in Holocene temperatures. Permafrost ice wedge δ¹⁸O, a clear winter season proxy, shows a continuous warming trend in the last 7 kyr BP (Meyer et al., 2015b) that supports the hypothesis that seasonality in proxy-based records is one key variable to be considered. Pollen-based reconstructions for the Holocene show a clear temperature decrease since an early to mid-Holocene temperature optimum (Fig. 8d; Herzschuh et al. 2021, 2022, 2023a, b) that is valid not only for the summer season but (though less pronounced) also for annual reconstructions. Proxy-based reconstructions of Kaufman et al. (2020) also suggest late Holocene cooling. The δ¹⁸O_BSi compilation presented here for lacustrine environments is based on diatoms whose bloom is mostly attributed to the late spring–early summer season.

Decreasing summer temperatures throughout the Holocene would manifest in decreasing T_lake, which would in turn lead to increasing δ¹⁸O_BSi (Fig. 1). The observed decrease in δ¹⁸O_BSi, however, points towards lower δ¹⁸O_lake, which would be in line with decreasing δ¹⁸O_prec due to decreasing summer T_air (Fig. 1). This would imply either a prevalence of summer precipitation or lake basins with sub-annual t_res. However, the trend is observed for lake basins with a wide range of t_res, suggesting millennial-scale changes in summer T_air to be the main driver of the observed δ¹⁸O_BSi signal.

It has to be stressed that this influence of T_air may act both directly and indirectly: directly by means of the temperature-dependent fractionation during precipitation formation and indirectly via the impact of temperature within the hydrological cycle and δ¹⁸O of water compartments. This includes factors such as moisture origin, precipitation intermittency, and atmospheric circulation patterns, which have been described as key drivers in numerous case studies. This underlines the importance of scale when assessing δ¹⁸O_BSi records, both temporally and spatially. On a millennial and hemispheric scale, T_air can be identified as one main driver of δ¹⁸O_BSi records even though the signal formation is generally complex and challenging to interpret with individual records. A T_air influence for certain records on the millennial scale is therefore not necessarily a contradiction to the findings of the original publications that may have attributed the record's signal to other factors (e.g. hydrological processes) on shorter timescales.

4.4 Frontiers and challenges (records older than the Holocene)

Further back in time than the Holocene, coverage with δ¹⁸O_BSi records is generally limited. As MIS 1 (N = 37) and MIS 2 (N = 18) have the highest numbers of records, they offer the best possibility of comparing glacial and interglacial records. Records beyond MIS 2 are very sparse and do not allow for filtering and generating common trends. Here, the interpretation has to rely on comparison of individual datasets. There are few (N = 16) records that cross the MIS 1–MIS 2 boundary, and these records cover neither MIS 1 nor MIS 2 entirely. A quantitative comparison of MIS 1 and MIS 2 is therefore difficult. Again, caution has to be applied due to the fact that these records stem from lakes with different hydrological and climatic settings (i.e. maritime Alaskan sites vs. continental Siberian sites).

The most promising approach for investigating MIS 1 and MIS 2 is extending NH MIS 1 datasets to the past (Fig. 9). Five NH δ¹⁸O_BSi records (#3, 8, 15, 24, 25) extend back from MIS 1 into MIS 2 and offer insight into the hydroclimate history of Siberia during the last glacial and the deglaciation. Four of these records (#3, 8, 15, 24) are continuations of records included in the NH stack (see Fig. 6f), and all of them stem from Asia. Most records display a maximum near the MIS 1–MIS 2 boundary and a decrease in δ¹⁸O_BSi with increasing age in MIS 2, suggesting decreasing T_air. This is the case for all records except #3, which displays a relative δ¹⁸O_BSi maximum at ca. 16 kyr BP. This maximum may also be caused by dry conditions outperforming the effect of lower T_air. Only records #3, #8, and #15 extend beyond 13 kyr BP. Both #3 and #15 reach their absolute minima at 23 and 24 kyr BP, respectively. Record #8 does not show such a clear pattern and shows a relative maximum at 25 kyr BP instead. This is likely due to the complex hydrological setting (alternating between open and closed at present) which may have changed over these timescales. Generally, most records suggest a tendentially lower δ¹⁸O_BSi in MIS 2 when compared to MIS 1.

Figure 9 Comparison of MIS 1 and MIS 2 data. Continuation of Holocene records from the NH stack (see Fig. 6f) into MIS 2 for colour-coded records (#3, 8, 15, 24).

This is notable as glacial and interglacial periods are characterized by different environments, atmospheric circulation patterns, and likely hydrological settings (e.g. formation or closure of outflows from lakes, lake level fluctuations). Potentially, lower δ¹⁸O_BSi values would be consistent with either lower T_air or more humid conditions. In the case of MIS2, generally associated with cold and dry conditions, a lower T_air is the more plausible scenario. However, it must be stressed that the lack of records covering the complete MIS 2 complicates a robust statistical comparison between MIS 1 and 2. Calculated offsets of MIS 1 and MIS 2 data points may show both higher and lower δ¹⁸O_BSi values in MIS 2 compared to the Holocene. Again, a MIS 2 climate colder and drier than Holocene (MIS 1) may produce either lower or higher δ¹⁸O_BSi values due to opposing effects of T_air and evaporation. A geographical pattern of either of these effects prevailing is not obvious, which suggests that rather than the climatic background, the individual hydrological settings of the lakes play a prominent role in determining the δ¹⁸O_BSi signal. This supports the approach of using lakes with both similar latitudinal and hydrological characteristics for meaningful inter-site comparison.

Beyond MIS 2, δ¹⁸O_BSi records become even scarcer, as do lake sediment records in general. Additionally, lake sediments covering these time periods often lack sufficient diatoms, especially in cold stages, i.e. Lake Baikal during MIS 4 (Mackay et al., 2008; Mackay et al., 2011; Mackay et al., 2013).

Lake El'gygytgyn is a peculiar example of a continuous sedimentation history, with δ¹⁸O_BSi showing glacial–interglacial cycles at least back to MIS 9 (Chapligin et al., 2012b), which have been attributed by the authors to T_air changes. In addition to glacial–interglacial cycles, it is also possible to compare δ¹⁸O_BSi of interglacials, i.e. MIS 1, MIS 5, MIS 11, when diatoms during cold stages are absent. In the sedimentary records of Lake El'gygytgyn, Chapligin et al. (2012c) investigated the δ¹⁸O_BSi differences between the warm stages MIS 1, MIS 5, and MIS 11 and found MIS 11 as warmest interglacial. This interpretation has been supported by pollen-based T_air reconstructions by Melles et al. (2012) and underlines the applicability of the δ¹⁸O_BSi proxy on glacial–interglacial timescales. Long-term diagenetic effects have shown to be of little influence on δ¹⁸O_BSi, at least for the last 250 kyr (Chapligin et al., 2012b). Studies on Lake Baikal have also investigated past interglacials and have addressed changes in precipitation intermittency and cooling events during MIS 11 (Mackay et al., 2008). During MIS 5, millennial-scale variability is suggested to have been more stable than during MIS 1 (Mackay et al., 2013).

Older δ¹⁸O_BSi records do exist and often rely on paleo-lakes, i.e. Ribains Maar, Baringo–Bogoria basin, Makgadikgadi, and Dethlingen (Shemesh et al., 2001b; Koutsodendris et al., 2012; Wilson et al., 2014b; Schmidt et al., 2017). These offer the possibility of extending the picture past extant lakes. However, the absence of modern hydrological observations makes their interpretation and comparison to other, better-constrained records more challenging. Further research using δ¹⁸O_BSi is therefore needed in order to complement the picture and provide insight into past climate and environmental conditions in continental regions, particularly valuable for high-latitude and high-altitude areas that are poorly covered by other proxy data. We recommend studies with uniform spatial and temporal coverage, i.e. hydrologically open lakes with a long, continuous sedimentation history, as our synthesis indicates these to be most promising for generating comparable binned time series to extend climate reconstructions further into the past.

5 Conclusions

In this study, we have identified and compiled all existing lacustrine δ¹⁸O_BSi records published to date and synthesized them into the first global δ¹⁸O_BSi data compilation. We have identified 53 δ¹⁸O_BSi records derived from 71 publications. These records stem from across the entire globe with their geographical distribution focusing on high-altitude and high-latitude lacustrine environments. Diatoms bear the advantage of being available in dilute, non-alkaline lakes common in these environments. Regional clusters of records (e.g. northern Eurasia, eastern Africa) indicative of the research foci of the individual research groups employing the δ¹⁸O_BSi proxy allow for generating regional subsets of data. Temporal coverage stretches from sub-recent timescales to the Pliocene, which highlights the applicability of δ¹⁸O_BSi on different timescales. Best coverage is available in the Holocene (MIS 1), which is linked to the age of lakes and hence the availability of lake sediments, especially at high latitudes. Moreover, biogenic silica (here entirely from diatom-based samples) is most abundant during warm periods (interstadials and interglacials) compared to colder periods. The interpretability of δ¹⁸O_BSi records relies on regional setting and hydrological constraints of individual lakes and catchments and is preferably supported by isotope measurements of modern lake water. Most δ¹⁸O_BSi records stem from open lakes (N = 41), suggesting for these lakes a negligible influence of lake evaporation. In contrast, closed lakes (N = 12) and paleo-lakes (N = 9) have been investigated less. It has to be noted that a lake's ontogeny and hydrology may change throughout time, and thus hydrological changes constitute a valid approach for interpreting these δ¹⁸O_BSi records.

Spatial resolution of the δ¹⁸O_BSi records is determined by the size of the lake and its corresponding catchment with lake water effectively integrating the input signal. While this signal may integrate large areas such as at Lake Baikal, most catchments (N = 18) are < 100 km², suggesting a locally confined signal for these lakes. Regarding temporal resolution, most records feature sampling resolutions which by far exceed t_res, suggesting sampling resolution to be the decisive factor in determining a record's signal. However, in case of t_res > 100 yr, lakes may be subject to increased lake evaporation.

Accounting for offsets and different temporal resolution of records and filtering for similar hydrological settings (hydrologically open lakes with t_res < 100 yr), we find a common pattern throughout the Common Era at the centennial scale, which we attribute to changing hydroclimate conditions. Changes are similar between Eurasia and North America, but they still differ in the timing of δ¹⁸O_BSi maxima. Generally, the combined Eurasian δ¹⁸O_BSi record seems to include major climate episodes during this period, including Roman Climate Optimum, Migration Period Pessimum, Medieval Climatic Anomaly (MCA), and Little Ice Age. An effect of recent warming is not visible in the Common Era data.

For the entire Holocene in Eurasia and the NH (for hydrologically open lakes with t_res < 100 yr), we find a common decreasing δ¹⁸O_BSi trend of ca. −0.2 ‰ kyr⁻¹ on the millennial scale which roughly follows summer insolation. The decreasing δ¹⁸O_BSi trend is accelerated in the second half of the Holocene, thus in Neoglacial times. Hence, this summer-based proxy contributes to the understanding of the Holocene climate history and the Holocene conundrum discussion initiated by Liu et al. (2014). The timing of the absolute δ¹⁸O_BSi maxima of individual records in Eurasia differs and shows an east–west gradient with eastern Eurasian records featuring earlier Holocene maxima compared to western Eurasian records. Holocene minima occur within the last 2 kyr for all Eurasian records. North American records diverge from this Holocene pattern with later maxima and earlier minima. This behaviour is likely linked to atmospheric circulation patterns as also described by the authors of the individual records.

Extending the millennial-scale trend of Eurasian records into MIS 2 shows generally lower δ¹⁸O_BSi values in MIS 2, suggesting lower T_air (and not $P/E$) in glacial times, also visible in most other records that act on glacial–interglacial timescales. The applicability of the δ¹⁸O_BSi proxy beyond MIS 2 is generally constrained by the availability of lake sediments and the abundance of diatoms in the respective sediments. Lake El'gygytgyn, as a prime example of hydrological continuity, displays glacial–interglacial cycles in the δ¹⁸O_BSi record, supporting T_air being a prominent driver of δ¹⁸O_BSi on millennial timescales. In glacial periods, diatoms abundances are low, and this limitation sometimes does not allow for δ¹⁸O_BSi analysis. However, comparison between interglacials where diatoms are generally more abundant is feasible, i.e. at Lake Baikal and Lake El'gygytgyn. In summary, we demonstrate the applicability and inter-comparability of combined δ¹⁸O_BSi records into a first lacustrine δ¹⁸O_BSi compilation allowing for paleoclimate reconstructions and accounting for regional offsets, different temporal resolutions, and hydrological backgrounds. δ¹⁸O_BSi records compiled in this study are an important tool for reconstructing paleoclimate across the globe and across a variety of timescales. In NH extratropical regions, their T_air-driven quantitative signal makes them especially useful in conjunction with paleoclimate models. Future research would be most valuable in complementing the existing data with longer records (covering MIS 2 and beyond) and records from underrepresented regions such as the Southern Hemisphere.

Appendix A

Figure A1 Geographical distribution of published δ¹⁸O_BSi records and year of publication (N = 64).

Figure A2 Overview of the dating methods used for age model creation of published δ¹⁸O_BSi time series (N = 64). Note that a single record can rely on multiple methods.

Figure A3 Overview of isotope hydrological background information considered in publications presenting δ¹⁸O_BSi records (N = 64).

Figure A4 Hydrological settings of the lakes corresponding to the δ¹⁸O_BSi records compiled in this study, as indicated in the original publications (N = 64). Owing to a lack of present constraints, paleolakes are listed as a separate category.

Figure A5 Holocene δ¹⁸O_BSi dataset for the African continent. (a) Original records, (b) data binned to 1 kyr intervals only showing records covering at least seven bins, (c) binned records with the mean of individual records removed, (d) with data including hydrological background information.

Table A1 Overview of the records and corresponding publications identified by this study (single meaning single core; composite meaning several cores; supplemented meaning information from other publication added, e.g. the age model).

Table A2 Overview of the records compiled and analysed in this study. Lake basin parameters shown are taken from the HydroLakes database (Messager et al., 2016) for consistency. In case of major discrepancies between the HydroLakes database and the values provided by individual case studies, values of case are given as well and marked with an asterisk.