Shireet Naiman Nuur is a high-altitude lake in the central Mongolian Khangai Mountains (46°31′55.04″N, 101°49′16.23″E; Figure 1). The lake is located at 2,429 m a.s.l. in the upper Orkhon Valley, where nomadic pastoralism has a long history (Figures 1, 2A). The lake’s small hydrological catchment (∼32.5 km²) with steep slopes reaches altitudes up to ∼3,100 m a.s.l. Shireet Naiman Nuur has a length of ∼3 km, a width of ∼3.4 km and a maximum water depth of ∼22.4 m. Currently, the lake has no outlet, i.e., it is endorheic, and the lake level is ∼5 m below a potential overspill into the upper Orkhon Valley. The regional climate in the Khangai Mountains and at Shireet Naiman Nuur is characterized by long and cold winters dominated by the Siberian High, whereas 75% of the moisture is supplied by the mid-latitude Westerlies during the short summer from June to August (Tian et al., 2013; Chen et al., 2019). Today’s mean annual temperature and precipitation at the nearest climate station Arvaikheer (CDC-ID 44388) is 2.0 ± 0.9°C and 256.4 ± 75.2 mm, respectively (period from 2009-2017; DWD Climate Center, 2021; Figure 2B). Because Shireet Naiman Nuur is located ∼600 m higher than the nearest climate station at Arvaikheer, an even shorter growing season with lower mean annual temperatures can be expected. Consequently, the extreme environmental conditions at the lake lead to an ice cover prevailing 7 months per year from November until May (2018-2020; Planet Team, 2017). The modern vegetation around the lake is dominated by grasses and steppe vegetation on the south-facing slopes, while the north-facing slopes are vegetated by relatively dense Larix sibirica forest because of greater moisture availability (Figure 2C).

We retrieved a 178 cm long gravity core in summer 2019 from the deepest part of Shireet Naiman Nuur using an Uwitec gravity corer. The sediment core was split and photographed at the Friedrich Schiller University Jena. One half of the core was analyzed for elemental composition in 2 mm resolution using an X-ray fluorescence (XRF) scanner (Avaatech) at Eawag, Switzerland. The XRF scanner was sourced by an Oxford 100 W X-Ray tube with a rhodium anode and equipped with a Canberra X - PIPS/DSA 1000 (MCA) detector. The split core surface was carefully flattened for uniform scanning and covered with 4 µm Ultralene foil for avoiding contamination. Low elemental groups were measured at 10 kV (1500 A, no filter, 15 s exposure) and mid energy groups were measured at 30 kV (2000 A, Pd thin filter, 40 s exposure). Elemental data from XRF-scanning is presented as log ratios to eliminate sediment matrix effects, i.e., due to water content, surface roughness and grain size variations (Weltje and Tjallingii, 2008). After XRF scanning, a u-channel was taken from the identical core half for paleomagnetic analyses. The other half of the sediment core was subsampled at 1 cm intervals, freeze-dried and homogenized.

The chronology of the sediment core is based on ¹⁴C-dating of nine bulk organic carbon and 23 terrestrial macrofossil samples. For collecting terrestrial macrofossils, 1 cm slices of the half of the sediment core where the u-channel was taken, were wet sieved over a 200 µm sieve and terrestrial macrofossils were subsequently collected and identified. Prior to ¹⁴C-analyses, macrofossils and bulk organic carbon samples were treated with 1N HCl at 60°C for 8 h to remove carbonates. ¹⁴C-dating was carried out at the Laboratory for the Analysis of Radiocarbon with AMS at the University of Bern (Szidat et al., 2014), using a MIni CArbon DAting System (MICADAS) AMS coupled online to an Elementar Analyzer (Ruff et al., 2010; Wacker et al., 2010). ¹⁴C-results are reported as fraction modern (F¹⁴C) and were corrected for cross and constant contamination after Salazar et al. (2015). 10 tin boats were measured for constant contamination, which gave 1.5 µg C for a single boat with F¹⁴C values of 0.498. ¹⁴C-ages were calibrated to calendar ages using the IntCal20 calibration curve (Reimer et al., 2020). ¹⁴C-ages were then used to establish an age-depth model by Bayesian age modeling using the Bacon 2.3.4 package in R (Blaauw and Christen, 2011). To account for possible reservoir effects of the bulk organic carbon ¹⁴C-ages, we used the uppermost bulk organic carbon sample from 0–1 cm sediment depth, which gave a conventional ¹⁴C-age of 632 ± 101. This reservoir age was subtracted from all bulk organic carbon ¹⁴C-ages during Bayesian age modeling. One bulk organic carbon sample at 80 cm sediment depth and two terrestrial macrofossils at 29 and 79 cm sediment depth gave too old ¹⁴C-ages and were therefore excluded from Bayesian age modeling in a second iteration.

For paleomagnetic secular variations, natural remanent magnetization (NRM) was measured on the u-channel using stepwise (11 steps) alternating field (AF) demagnetization with peak AFs of 0, 5, 10, 15, 20, 30, 40, 50, 65, 80 and 100 mT. NRM was measured continuously at 1 cm intervals using a 2G Enterprise DC-4K liquid helium free magnetometer at GFZ Potsdam. The u-channel was split into two halves during sampling and paleomagnetic measurements on the u-channels were performed in 1 cm steps extending 10 cm over both ends of the u-channel profiles with respect to the absolute depth scale. These data sets, on basis of the measured cartesian components (X,Y,Z), were then superimposed to compensate for the end-effect at the edges of the U-channels, thus yielding a continuous paleomagnetic record (Nowaczyk et al., 2020). Inclination and declination of the characteristic remanent magnetization (ChRM: in most cases 15–80 mT) and maximum angular deviation (MAD) were determined by principal component analyses (Kirschvink, 1980) and the median destructive field was also calculated. Declination data are relative since the azimuth could not be controlled during coring. Due to the width of the response function of the SQUID sensors, some smoothing occurs in the data. In order to eliminate edge effects, the upper- and lowermost 5 cm of the sequence are plotted in gray (Figure 3B).

Analyses for carbon and nitrogen contents, and the stable carbon isotopes were carried out with an Elementar Analyzer (vario EL cube) coupled to an isotope ratio mass spectrometer (Isoprime precision) in 1 cm resolution at the Friedrich Schiller University Jena. For total carbon (TC) and total nitrogen (N), ∼10 mg of untreated sediment samples was weighed into tin boats and measured. For total organic carbon (TOC) and the carbon isotopic composition (δ¹³C_TOC), samples were treated with 1N HCl at 60°C for 8 h to remove carbonates. The samples were subsequently washed to pH neutrality with ultrapure water, and ∼10 mg of sediment was weighted into tin boats for measurement. The analytical precision of the δ¹³C analyses was checked against certified standards (L-Prolin, EDTA and USGS65), and gave an analytical error < 0.1‰. δ¹³C is given in its delta notation against the Vienna Pee Dee Belemnite (VPDB).

The molar TOC/N ratio was calculated based on their respective molecular weights as follows: TOC / N [ molar ] = ( TOC [ % ] ∗ 12 . 0107 [ g∗mol - 1 ] ) / ( N [ % ] ∗ 14 . 0067 [ g∗mol - 1 ] ) (1)

Pearson’s r-values were calculated to detect correlations between the analyzed geochemical proxies and a student’s t-test was used to check for significance (α < 0.05). Statistical analyses were carried out with the software OriginPro 2022.

The 178 cm long sediment core from Shireet Naiman Nuur comprises mostly layered sediments of fine- to medium-grained silty material (see Supplementary Information SI 2 for grain size analyses and results) and shows three lithological units based on differences in sediment color and geochemical properties. In Unit A (from 178 to 118 cm), two pale layers at 167 and 152 cm sediment depth interrupt the brownish and well laminated sediments. The brownish layered sediments in Unit B (from 118 to 60 cm) are interrupted by a double pale layer at 114 cm sediment depth, and blackish sediments at 80 cm sediment depth. In Unit C (from 60 to 0 cm), another pale layer occurs at 59 cm sediment depth and well laminated brownish sediments contain intercalated blackish layered sediments (Figure 3A).

Bayesian age-depth modeling of ¹⁴C-dated bulk organic carbon samples and terrestrial macrofossils (see Supplementary Information SI 3 for detailed ¹⁴C-dating results) gave a modeled basal median ¹⁴C-age of 7.4 ± 0.3 cal ka BP. All ¹⁴C-ages are in stratigraphic order, except for a terrestrial macrofossil sample at 29 cm, which shows a large offset to the other ages, as well as a bulk organic carbon sample at 80 cm and another terrestrial macrofossil sample at 79 cm, which are slightly too old for their stratigraphic position but still within the modelled 2σ error (Figure 3A). The reservoir corrected bulk organic carbon ¹⁴C-ages correspond very well with the ¹⁴C-ages of the terrestrial macrofossils.

NRM intensities range from 2.4 to 28.8 Am⁻¹ and MDF values from 27.9 to 43.3 mT. Inclination and declination show distinct changes (Figure 3B). Inclination ranges from 52 to 65°N and only intersects the one expected, based on a geocentric axial dipole model (GAD = 64.64°), in the lowermost part of the record. Such a flattened inclination compared to the GAD has been often noticed in sediments and has been attributed to compaction processes (Henkel et al., 2016). Declination stretches over 32°, which in general is typical for paleomagnetic secular variations (Haberzettl et al., 2015; Haberzettl et al., 2021). Well preserved and stable single-component magnetization is indicated by low MAD values (Stoner and St-Onge, 2007) that are generally below 4.1°.

TOC and N contents range from 4.5 to 10.6% and 0.6–1.4%, respectively, with generally higher values in Unit A (i.e., from 178 to 118 cm), a decreasing trend in Unit B (i.e., from 118 to 60 cm) and lower values in Unit C (i.e., from 60 to 0 cm; Figure 4). The molar TOC/N ratio shows values between 5.4 and 7.9, covering only a very small range. The pattern of TOC and N is resembled by log (Ca/Ti) and log (Si/Ti) values that are likewise higher in Unit A, show a decreasing trend in Unit B and lower values in Unit C. Abrupt decreases in TOC, N, log (Ca/Ti) and log (Si/Ti) are notable at 167, 152, 114 and 59 cm sediment depth, corresponding to the four macroscopically detected pale layers (Figure 4). δ¹³C_TOC values range from -31.4 to -27.9‰ and are more negative in Unit A, variable in Unit B and more positive in Unit C (Figure 4).

Correlations of the geochemical proxies are all significant (Table 1). TOC, N, log (Ca/Ti) and log (Si/Ti) are positively correlated with each other, but negatively correlated to δ¹³C_TOC.

Extensive ¹⁴C-dating of the finely layered silty sediments from Shireet Naiman Nuur allows us to establish a robust and precise chronology for the past 7.4 ± 0.3 cal ka BP (Figure 3A). While almost all ¹⁴C-ages are in stratigraphic order, the terrestrial macrofossil sample at 29 cm sediment depth is pre-aged and has an age-offset of ∼1,500 years. Although the median ¹⁴C-age of the bulk organic carbon sample and the terrestrial macrofossil sample at 80 and 79 cm sediment depth are still within the modelled 2σ error, they have a slight age-offset of ∼250 years and are likewise slightly pre-aged (Figure 3A). Those pre-aged samples probably represent reworked and mobilized organic material that has been stored in the catchment soils for a certain time before final erosion (Haberzettl et al., 2013; Gierga et al., 2016; Bliedtner et al., 2020). They therefore overestimate their “true” burial age and we excluded them from age-depth modeling. However, all remaining ¹⁴C-ages are stratigraphically very consistent and indicate that residence times of organic material in the catchment soils and transfer times through the catchment are very short and negligible. Moreover, the temporal resolution of our core derived from Bayesian age modeling is 16 (i.e., in the upper part) to 72 years (i.e., in the lower part) per centimeter sediment, which correspond to a decadal to centennial temporal resolution.

The ¹⁴C-based chronology from Shireet Naiman Nuur is further confirmed by paleomagnetic secular variation stratigraphy. Observed variations in inclination correspond very well with outputs of various spherical harmonic geomagnetic field models. A comparison to the pfm9k model (Nilsson et al., 2014) is shown in Figure 5, which showed the best resemblance. Slight offsets are within the ¹⁴C-dating uncertainties of the Shireet Naiman Nuur sediments. Compared to the closest regional paleomagnetic secular variation record from Lake Baikal (Peck et al., 1996), the Baikal inclination pattern shows similar features to Shireet Naiman Nuur and the pfm9k inclination pattern (Figure 5A), but less obviously and sometimes with an offset in timing. Nevertheless, the overall trend of the inclination pattern at Lake Baikal follows that of our record and the pfm9k model during the past 7.4 ka. Discrepancies in the inclination pattern of Lake Baikal compared to Shireet Naiman Nuur and the pfm9k model might be explained by small-scale regional differences in the drift of the non-dipole field (Merrill et al., 1996). However, they could also be explained by possible uncertainties in the Lake Baikal chronology, that is, mostly based on ¹⁴C-dated bulk organic carbon (Colman et al., 1996; Peck et al., 1996). Although the ¹⁴C-ages are reservoir corrected by dated surficial sediments, changes in the discharge of the Selenga River into Lake Baikal could have strongly influenced sedimentation rates and transfer times of organic material deposited in the analyzed sediment cores located in front of the Selenga delta. Comparing the declination pattern recorded in the Shireet Naiman Nuur and Lake Baikal sediments over the past 7.4 ka, both records resemble very well the declination pattern of the pfm9k model (Peck et al., 1996; Nilsson et al., 2014). Especially since ∼4 ka, three distinct declination lows and highs are recorded, although H1 and L1 appear a bit earlier in the Shireet Naiman Nuur and Lake Baikal record compared to the pfm9k model (Figure 5B).

Overall, inclination and declination patterns recorded in the Shireet Naiman Nuur sediments are in very good agreement with the pfm9k model of Nilsson et al. (2014). Therefore, paleomagnetic secular variations validate and confirm our extensively dated ¹⁴C-chronology, providing a very robust age control for the 7.4 ka lake sediments from Shireet Naiman Nuur. Moreover, our robustly dated paleomagnetic secular variation record can serve as a valuable regional master record for using paleomagnetic secular variation stratigraphy in semi-arid Mongolia as a powerful tool to evaluate and refine ¹⁴C-chronologies. This is especially important since terrestrial macrofossils are often absent in semi-arid lake sediments and bulk organic carbon can have large residence and transfer times in and through the catchment.

Changes in TOC in the Shireet Naiman Nuur sediments are interpreted in terms of lake primary productivity. The higher TOC values in Unit A followed by the trend to lower values in Unit B and the lower values in Unit C indicate a trend from increased to low primary productivity in the Shireet Naiman Nuur sediments (Figure 4). TOC is a well-known indicator for primary productivity in lacustrine sediments if organic material is mainly derived from aquatic biogenic producers (Haberzettl et al., 2005; Kasper et al., 2015; Zhu et al., 2015), which is suggested for the Shireet Naiman Nuur sediments by low molar TOC/N ratios < 8 (Meyers, 2003). TOC is additionally resembled by N, log (Ca/Ti) and log (Si/Ti), and their significant correlation with TOC (r = 0.91, 0.65 and 0.51, respectively; α < 0.05; Table 1) suggesting that they are also indicative for changes in lake primary productivity in the Shireet Naiman Nuur sediments. We acknowledge that Ca and Si could also originate from geogenic sources from the catchment or productivity independent precipitation processes. However, due to the significant correlations to TOC and N, log (Ca/Ti) and log (Si/Ti) ratios are similarly interpreted. Changes in primary productivity in high-altitude and high-latitude lakes in turn are mainly controlled by growing season temperatures because the duration of ice-cover and the growing season of aquatic producers is strongly determined by air temperatures (Willemse and Törnqvist, 1999; Mischke et al., 2010). Recently, such a relationship was also found for the past 4.2 ka in the high-altitude lake sediments of Khar Nuur which is located in the Mongolian Altai and has a comparable environmental setting (Bliedtner et al., 2021). According to this interpretation, lake primary productivity and growing season temperatures at Shireet Naiman Nuur were high in Unit A from 7.4 ± 0.3 to 4.3 ± 0.2 cal ka BP, decreasing in Unit B from 4.3 ± 0.2 to 1.9 ± 0.2 cal ka BP and lowest in Unit C from 1.9 ± 0.2 cal ka BP until present day (Figures 4, 6). The abrupt decreases of lake primary productivity in the pale layers at 167 cm (6.9 ± 0.3 cal ka BP), 152 cm (6.3 ± 0.2 cal ka BP), 114 cm (4.2 ± 0.2 cal ka BP) and 59 cm (1.8 ± 0.2 cal ka BP) sediment depth probably document increased minerogenic input from erosive events, possibly induced by rapid climate changes (Figures 4, 6).

Lower growing season temperatures and longer ice cover periods in Unit C (i.e., from 1.9 ± 0.2 cal ka BP until present day) are supported by more positive δ¹³C_TOC values (Figure 4). Since organic material in the Shireet Naiman Nuur sediments is mainly derived from aquatic biogenic production (see low molar TOC/N ratios), the δ¹³C_TOC signal mainly indicates the dissolved inorganic carbon (DIC) sources aquatic producers use for biosynthesis. δ¹³C_TOC values are more negative if aquatic producers use dissolved CO₂ from the water column, that is, in equilibrium with the atmosphere (δ¹³C = ∼−7‰). In contrast, δ¹³C_TOC values become more positive when dissolved CO₂ is limited and exhausted, and aquatic producers begin to use dissolved bicarbonates (HCO₃ ⁻) (δ¹³C = ∼1‰) (Meyers, 2003). In most cases enhanced primary productivity results in the consumption of dissolved CO₂ and the incorporation of dissolved HCO₃ ⁻ (Meyers, 2003; Leng and Marshall, 2004; Lamb et al., 2006), but in high-altitude and high-latitude lakes also ice cover can determine the availability of dissolved CO₂ (Striegl et al., 2001; Neumann et al., 2004). Studies of high-latitude lakes showed that a shortening of the ice cover period results in longer mixing of the lake during spring with a slower rate until summer stratification sets in, and vice versa if the ice cover period was longer (Couture et al., 2015; Pilla and Williamson, 2021). Therefore, we suggest at Shireet Naiman Nuur, where ice coverage is in the order of 7 months, that more dissolved CO₂ is available during shorter ice periods with longer mixing, resulting in more negative δ¹³C_TOC values, i.e., from 7.4 ± 0.3 to 4.3 ± 0.2 (Figure 6). In contrast, dissolved CO₂ becomes exhausted during longer ice cover periods with shorter mixing, resulting in beginning incorporation of dissolved HCO₃ ⁻ and more positive δ¹³C_TOC values, i.e., from 1.9 ± 0.2 cal ka BP until present day. Variable δ¹³C_TOC values from 4.3 ± 0.2 to 1.9 ± 0.2 cal ka BP indicate changing dissolved DIC sources due to changes in ice cover periods (Figure 6). Lower growing season temperatures and longer ice cover periods are further supported by the significant negative correlation of the lake primary productivity indices TOC and log (Ca/Ti) with δ¹³C_TOC (r = -0.60 and -0.63, respectively; α < 0.05; Table 1).

Therefore, increased primary productivity and higher growing season temperatures occurred at Shireet Naiman Nuur during the Middle Holocene from 7.4 ± 0.3 to 4.3 ± 0.2 cal ka BP. Lower growing season temperatures and longer ice cover periods occurred during the Late Holocene from 4.3 ± 0.2 cal ka BP until present day (Figure 6).

The overall pattern of higher growing season temperatures at Shireet Naiman Nuur during the Middle Holocene seem to be driven by higher summer insolation (Figure 6). Strong summer solar heating, Arctic temperature amplification and sea ice loss was reported for mid- and high-latitudes of the Northern Hemisphere for the Middle Holocene, resulting in the “Holocene thermal maximum” and a pronounced warming (Wanner et al., 2008; Renssen et al., 2009; Park et al., 2019). Especially for mid-latitude Central Asia, increased summer temperatures between ∼8 and 4 ka BP were recently reported to be insolation-driven by a compiled temperature index of Lan et al. (2021). With the onset of neoglacial cooling ∼4 ka BP due to reducing summer insolation, colder temperatures prevailed during the Late Holocene (Zhao et al., 2017; Lan et al., 2021), which is also evident at Shireet Naiman Nuur during the past 4.3 ± 0.2 cal ka BP (Figure 6). Although growing season temperatures overall started to decrease since then, an abrupt decrease in growing season temperatures at 4.2 ± 0.2 cal ka BP is notable and exceptional in our record. This abrupt cooling likely indicating an event of rapid climate change at Shireet Naiman Nuur and coincides with the 4.2 ka event, a period of cold conditions probably related to enhanced ice-rafted debris in the northern North Atlantic (Bond et al., 2001). Unfortunately, it was not yet possible to attribute the other pronounced minima in our productivity indicators at 6.9 ± 0.3 and 6.3 ± 0.2 cal ka BP to more widely known climate anomalies.

While summer insolation generally decreased during the Late Holocene, variations in growing season temperatures at Shireet Naiman Nuur during that time might be modulated by total solar irradiance (TSI) (Steinhilber et al., 2009). Although we have to note that the timing of growing season minima slightly differs from minima in TSI, they are still within our chronologies 2σ error. Moreover, the magnitude of changes in the productivity indices should rather be considered as relative changes because of site-specific and matrix effects, and therefore do not necessarily reflect the same magnitude of changes in TSI. Nevertheless, abrupt changes to lower growing season temperatures occur at Shireet Naiman Nuur during minima in TSI at ∼3.4, 2.8, 2.4 cal ka BP. A strong cooling starts ∼1.8 cal ka BP and lasts until ∼1 cal ka BP (Figure 6). This cooling phase includes the Dark Ages Cold Period (DACP) and/or Late Antique Little Ice Age (LALIA) from ∼1.6 to 1.2 ka BP (Büntgen et al., 2016; Helama et al., 2017). Compared to Mongolian lowland records (Yang et al., 2020; Rudaya et al., 2021; Struck et al., 2022), the following Medieval Climate Anomaly (MCA) does not appear that pronounced as a warm period in the Shireet Naiman Nuur sediments. However, colder conditions are evident during the Little Ice Age (LIA) and growing season temperatures especially decrease during the Maunder Minima (Figure 6).

The insolation-driven trend of warmer growing season conditions during the Middle Holocene to colder growing season conditions during the Late Holocene since 4.3 ± 0.2 cal ka BP is also reflected in several pollen reconstructions from central and northern Mongolia (Fowell et al., 2003; Prokopenko et al., 2007; Wang et al., 2011; Unkelbach et al., 2021). The TSI controlled variation in growing season temperatures at Shireet Naiman Nuur during the Late Holocene was previously described for another high-altitude site from the Mongolian Altai (Bliedtner et al., 2021). Although the MCA is not that pronounced in both records, they well reflect the Late Holocene climate anomalies that were also regionally reported from pollen reconstructions (Yang et al., 2020; Rudaya et al., 2021), a glycerol dialkyl glycerol tetraether based temperature reconstruction (Dugerdil et al., 2021) and an isotope based moisture reconstruction (Struck et al., 2022).