Hydroclimate in the Pamirs Was Driven by Changes in Precipitation‐Evaporation Seasonality Since the Last Glacial Period

The Central Asian Pamir Mountains (Pamirs) are a high‐altitude region sensitive to climatic change, with only few paleoclimatic records available. To examine the glacial‐interglacial hydrological changes in the region, we analyzed the geochemical parameters of a 31‐kyr record from Lake Karakul and performed a set of experiments with climate models to interpret the results. δD values of terrestrial biomarkers showed insolation‐driven trends reflecting major shifts of water vapor sources. For aquatic biomarkers, positive δD shifts driven by changes in precipitation seasonality were observed at ca. 31–30, 28–26, and 17–14 kyr BP. Multiproxy paleoecological data and modelling results suggest that increased water availability, induced by decreased summer evaporation, triggered higher lake levels during those episodes, possibly synchronous to northern hemispheric rapid climate events. We conclude that seasonal changes in precipitation‐evaporation balance significantly influenced the hydrological state of a large waterbody such as Lake Karakul, while annual precipitation amount and inflows remained fairly constant.


Introduction
The Pamirs are located at the intersection of atmospheric influences from the Asian monsoon system, the midlatitude westerly winds (westerlies), and the wintertime Siberian anticyclone (Aizen et al., 2001). They are characterized by a pronounced gradient of decreasing winter and spring precipitation from west to east (Miehe et al., 2001;Williams & Konovalov, 2008). Past changes in hydroclimate in this region were likely to experience shifts in the extent of those large-scale atmospheric circulation systems (Caves Rugenstein & Chamberlain, 2018;Zhang et al., 2017). There is a consensus that the Asian monsoon is driven by insolation on a precessional scale, with increasing intensity since the late glacial and a maximum in the early to mid-Holocene (Herzschuh, 2006;Wang et al., 2017). However, mechanisms and the local and regional expression Recently, a lacustrine archive with a robust chronology spanning a glacial-interglacial cycle was obtained from Lake Karakul . Inorganic geochemical data, sedimentological, and paleoecological proxies from this record showed large-scale shifts, possibly driven by insolation, over the past 28 kyr (Heinecke, Epp, et al., 2017;Heinecke et al., 2018). Some of those proxies were interpreted with respect to past lake level changes, but hydrological inferences were tentative.
To decipher past hydrological conditions in Lake Karakul and its catchment in more detail, we analyzed geochemical proxies including the hydrogen isotopic signatures of both terrestrial and aquatic biomarkers whose values are directly linked to the hydrological cycle. Results were placed in context to a set of experiments using two atmospheric models with embedded isotope modules CAM3-iso  and ECHAM5-wiso (Werner et al., 2011), as well as the climate model EC-Earth (Hazeleger et al., 2010;Hazeleger et al., 2012). With this coupled proxy-model approach and with consideration of prior ecological and geochemical evidence, we aim to (I) evaluate regional hydrological responses to large-scale climatic shifts over a glacial-interglacial cycle; and (II) decipher potential drivers behind inferred lake-level changes.

Materials and Methods
A detailed description of methodology is given in supporting information Text S1. Briefly, a composite agemodel was built for two parallel sediment cores (KK12-1 and -2) taken in 12-m water depth (39. 01760°N and 73.53275°E;Mischke et al., 2017). For organic geochemical analysis, samples from core KK12-2 were extracted with an organic solvent and split into fractions by using silica-gel columns. Concentrations, carbon, and hydrogen isotopic values of compounds in the aliphatic and ketone fraction were obtained by gas chromatography and isotope ratio mass spectrometry under application of standard measures for quality assurance. The aquatic to terrestrial ratio (P aq ) was calculated according to Ficken et al. (2000): P aq = (nC 23 + nC 25 )/(nC 23 + nC 25 + nC 29 + nC 31 ). Δ terr-aq was calculated as 1,000 × [(δD-terr + 1,000)/(δD-aq + 1,000) − 1] using the arithmetic mean between δD-C 23 and δD-C 25 as δD-aq and between δD-C 29 and δD-C 31 as δD-terr. Elemental data were derived from XRF-scanning and inorganic isotopic values obtained by methods described in Heinecke, Epp, et al. (2017),  and-if measured at core KK12-1 -were transferred to the new composite age-model.

Biomarker Abundance and Composition-Assessing Organic Matter Sources
The organic matter in most sections was dominated by aquatic sources (P aq value~0.8; Ficken et al., 2000;Figure 2i), evidenced by visible weakly decomposed macrophyte remains and high concentrations of nC 23and nC 25 -alkanes (approx. 1.5-12 μg/g d.w.) in the sediment extract, especially after ca. 7 kyr BP (compiled stack-plot of all measured parameters in the supporting information Figure S1). The relatively high δ 13 C signatures of those compounds (approx. −15‰) provide further support for a submerged aquatic macrophyte and macroalgae origin (Aichner et al., 2010) with Potamogetonaceae and Characeae as the predominant families (Heinecke, Epp, et al., 2017). Less abundant nC 29and nC 31 -alkanes (0-2.5 μg/g d.w.) with δ 13 C values of approximately −30‰ are considered to reflect the sparse alpine C 3 vegetation in the catchment today, possibly with minor contribution from biomarkers derived from species thriving in small patches of wetlands around the lake. Terrestrial compounds were present in lower concentrations (<0.5 μg/g) during most parts of the glacial as expected during glacial expansion and cooling, and a more persisting ice-cover on the lake. Three episodes revealed increased terrestrial plant wax concentrations during the late Marine Isotope Stage (MIS) 3 and early MIS 2 at ca. 31-30 and 28-26 kyr BP, and during the late glacial at 17-14 kyr BP ( Figure S1).  Taft et al., 2014). Seasonal climatic data (blue bar, mean monthly precipitation; red curve, mean monthly temperature) from Karakul station (yellow dot; 1934-1990) in comparison to Sarytash (ca. 60 km north;1934-2000 and Fedchenko glacier (ca. 65 km west;1936-1994Williams & Konovalov, 2008). Insert: A = Kesang Cave. B = Ton Cave. C = Darai Kalon loess section.  (Berger & Loutre, 1991). (b, c) Karakul (KK) long-chain leaf-wax δD values in comparison to Ti areas from XRF scanning and (f, g) δ 18 O data from Ton and Kesang Cave speleothems (Cheng et al., 2016) and (h) δD values from Darai Kalon loess section (Haggi et al., 2019). (e) Δ terr-aq : Offset between δD values of aquatic and terrestrial n-alkanes in comparison to (d) winter westerlies index (WWI; Li et al., 2013). Bottom: Proxies indicative for lake level and vegetation changes. (i) P aq : Proxy for aquatic macrophytes (Ficken et al., 2000). (j) Summer insolation. (k) Reconstruction of 15 kyr BP lake level highstand (Komatsu & Tsukamoto, 2015). (l) PCA axis 2, indicating lake level  based on inorganic proxies measured at core KK12-1, adjusted to the new composite age-model for both parallel cores KK12-1 and -2, as published in Mischke et al. (2017). (m, n) Inferred salinity changes based on Sr/Ca ratios and alkenone abundances. Those intervals also have a low aquatic to terrestrial ratio (P aq value~0.3, Figure 2i) and were further characterized by the occurrence of unsaturated compounds (nC 23 , nC 25 , and nC 27 alkenes; Figure 2q). While the exact origin of those molecules remains unclear, a source from phytoplankton such as eustigmatophytes and chlorophytes and possibly from insects is likely (Aichner et al., 2012;van Bree et al., 2014;Zhang et al., 2015). Occurrence of n-alkenes correlate with increased counts of Chironomid remains at 28-26 and 15-14 kyr BP (Heinecke et al., 2018). Conversely, a much broader range of compounds than observed in the Karakul sediments would be expected with insects as a significant precursor of those compounds (Chikaraishi et al., 2012;van Bree et al., 2014). Long-chain C 37 alkenones derived from phytoplankton (de Leeuw et al., 1979) were abundant between ca. 17-5 kyr BP (10-200 ng/g d.w; supporting information Figure S1). For the overall trend of organic matter input, bromine content derived from XRF scanning is used as an estimate ( Figure S1), because many organic molecules contain carbon bound bromine, which is released after decomposition (Mayer et al., 1981).

Biomarker Stable Isotopes as Recorders of Past Local Hydrology and Moisture Sources
Inorganic and organic geochemical proxies showed distinct glacial-interglacial trends over the past 31 kyrs (Figures 2 and S1). Since nC 29 and nC 31 alkanes were mainly derived from the sparse alpine meadows around the lake, their δD values reflect a precipitation signal recorded during the vegetation period (i.e., summer, with a possible contribution of a winter signal by early spring snow-melt) within the lower altitude catchment ( Figure 2c). In contrast, we interpret aquatic nC 23 and nC 25 alkane δD values to record a lake water isotopic signal that integrates over the whole year, including meltwater derived from higher altitudes ( Figure 2p). Because of the interpretations for δD-aq (winter + summer signal) and δD-terr (summer signal), we infer the offset (Δ) between them as a proxy for winter precipitation and thus being driven by the location and intensity of the westerly flow during winter (Figure 2e).
The higher δD-terr values during the glacial and lower values during the Holocene are in agreement with speleothem δ 18 O data from Ton Cave and Kesang Cave 8 (Cai et al., 2017;Cheng et al., 2016; Figure 1) and leaf wax δD values from a loess section in western Tajikistan (Haggi et al., 2019), highlighting the transregional relevance of this signal (Figures 2f, 2g, and 2h). Seasonal isotopic variability of precipitation in inner Asia is driven by temperature (Araguas-Araguas et al., 1998;Yao et al., 2013). However, the higher δD-terr values during the glacial are unlikely to reflect a higher temperature. Higher titanium counts during the glacial ( Figure 2b) could indicate drier conditions; however, the large isotopic amplitude within the δD-terr data (approx. 60‰) cannot be solely explained by further aridification of a lower lake basin already characterized by present day MAP of <100 mm. Hence, we suggest that the shift is driven by a large-scale change of water vapor sources between glacial and Holocene summers, while higher titanium abundance in this case indicates intensified winds rather than further aridification during colder conditions.
In agreement with this interpretation of titanium counts, the CAM3-iso model simulated higher winter windspeeds for LGM and H1 compared to the Holocene (Figures 3 and 4a). This is related to an LGM eastward shift of North Atlantic storm tracks that resulted from atmospheric responses to the expansion of ice-sheet and sea-ice (Kageyama, D'Andrea, et al., 1999;Laine et al., 2009). Wind fields further illustrate the dominance of the southern branch of westerlies during winter and spring and resulting southwesterly wind in the Pamirs over the studied time intervals ( Figure S2a). In contrast, in summers the region is more dominated by the northern branch with northern winds (in agreement with computed back trajectories by Yan et al., 2019), which is an additional explanation for the lower precipitation compared to winter and spring.
CAM3-iso also showed that most of the summer moisture in the study area is derived from continental recycling (supporting information Figure S3). Thus, the relevant driver behind isotopic variability during summers is the balance between (more enriched) local and continental moisture sources and those which have been transported farther (more depleted). The shift toward more depleted δD-precip values at the glacial-interglacial boundary could therefore be interpreted as a strengthening of either the monsoon (increased transport from Arabian Sea) or of the westerlies (increased vapor transport from North Atlantic; Figure S3).
Adding to this complexity is the observation that monsoonal strength and westerly air-masses are closely linked to each other. An intensified summer monsoon in East Asia is associated to a jet-stream located over  Table S1. Temperature and precipitation data in comparison to present day instrumental data (colored lines; Williams & Konovalov, 2008). δD data in comparison to data from the Online Isotope in Precipitation Calculator (OIPC; red line; Bowen & Revenaugh, 2003). RH = relative humidity; PME = precipitation minus evaporation.
the Pamirs (Schiemann et al., 2009), resulting in stronger summer westerlies. Based on a 150-kyr long simulation, Li et al. (2013) suggested that phases of increased northern hemispheric insolation and more intense Asian monsoon are potentially also accompanied by an increase of winter precipitation in Central Asia, depending on latitude. They proposed that low winter insolation and resulting high latitude cooling decrease the northern hemispheric temperature gradient, (low winter westerlies index WWI; Rossby, 1939), which in turn weakens the northern westerly streams but enforces the southern branch. It is difficult to assess the consequences of those findings with regard to precipitation amount changes in our study area, as it is located at the proposed latitudinal boundary of 40°N. However, our suggested indicator for westerlies intensity (Δ terr- Figure 4. Anomaly plots LGM versus MH for summer (JJA; left) and winter (DJF; right) for (a) windspeed (850-hPa sigma level, derived from CAM3-iso), (b) precipitation, and (c) precipitation minus evaporation (PME) both derived from EC-Earth. Study area indicated by arrow. aq ) is in phase with winter insolation (Berger & Loutre, 1991) and the WWI (Li et al., 2013), supporting this interpretation (Figures 2a, 2d, and 2e). Weak winter westerlies between ca. 17-7 kyr BP in our study area correspond to findings from the Tien Shan, where similar conditions have been reconstructed for the late glacial and early Holocene, based on grain size analysis of a loess section (Jia et al., 2018).
A late glacial and early Holocene intensification of the Asian monsoons is well established (Hou et al., 2017;Rao et al., 2016;Wang et al., 2005), but it is unclear if the Pamirs were significantly influenced by moisture derived from the Indian summer monsoon, even though this is suggested by the most negative δD-terr values in our Karakul record between ca. 12-7 kyr BP. The co-occurrence of strengthening of monsoons and westerlies in Central Asia makes it difficult to determine the precise source of moisture for those time intervals. The modelling data suggest that precipitation over the region is controlled by the westerlies all around the year and for the studied time periods (Figures S2a and S2b). On the contrary, based on the model outputs, it also becomes clear that the western Pamirs are affected in a significantly stronger way by the incoming summer north-westerlies than our study area in the eastern Pamirs ( Figure S2a), resulting in the observed E-W gradient of spring precipitation ( Figure S2b) and of summer PME (precipitation-evaporation balance; Figure S2c). This illustrates that the Karakul region is situated in a transitional area in terms of precipitation amount, water balance, and wind directions for all simulated time intervals.
We conclude that the changes of atmospheric circulation patterns during the late glacial to Holocene transition were characterized by a decreased winter westerly flow and a shift from local to distant vapor sources during summers. The stronger summer monsoon over southeast Asia and India during that time significantly influenced the regional hydroclimate in the Pamirs-either directly by the monsoon itself or more likely by influencing the pathways and intensity of the westerlies.

Lake Level Highstands During Cold Episodes Caused by Wetter Summers
Several independent proxies provide evidence that Lake Karakul reacted sensitively to changes of local hydrological and climatic conditions. Based on OSL dating of paleoshorelines, Komatsu and Tsukamoto (2015) identified several lake transgressions, the last starting before 19 kyr BP and forming a shoreline 35 m higher than present around 15 kyr BP (Figure 2k). In addition,  interpreted higher values on the second axis of a principal component analysis of geochemical parameters, described by high TIC content and Fe/Mn ratios, to reflect glacial and late glacial lake-level highstands (Figure 2l). In our sedimentary record, we further interpret an increasing Sr/Ca ratio after 14 kyr BP as indicator for increasing salinity in a shrinking lake (Dissard et al., 2010;Yan & Wünnemann, 2014;Figure 2m). The interval 14-10 kyr BP was also characterized by decreasing proportional abundance of tetra-unsaturated alkenones C 37:4 (Figure 2n), which was found to be related to higher salinities based on studies of surface sediment samples from lakes on the Tibetan Plateau (Liu et al., 2008;Liu et al., 2011) even though other factors such as temperature, nutrient availability, and changing producers might complicate this relationship.
Taken together, all these proxies suggest rising lake level in the late glacial (peaking around 15 kyr BP) and episodically also in the last glacial period (between ca. 31-30 and 28-26 kyr BP). These phases are accompanied by ecological changes within the lake (shift in dominant submerged plants at coring location; occurrence of biomarkers such as alkenones and C 23 , C 25 , and C 27 n-alkenes, giving evidence for a change in algal communities) as well as around the lake (change from desert to steppe vegetation; Heinecke et al., 2018). Unsaturated midchain compounds were discovered in many Tibetan Plateau lakes. Their precise source is still unknown, but analysis of surface sediments of Lake Donggi Cona revealed that they mainly occur in greater depths (Aichner et al., 2012). At present, Chara species occur rather in deeper parts of Lake Karakul, while Potamogetonaceae thrives down to a maximum depth of approximately 14-15 m (Heinecke, Epp, et al., 2017;. The massive occurrence of weakly decomposed Potamogetonaceae remains after ca. 7 kyr BP, and synchronous increase of nC 23 and nC 25 concentrations in sediments were therefore most likely induced by the combined effects of a falling lake level, change in lake water chemistry (disappearance of alkenones), and resulting changes in light penetration to the ground, leading to pronounced propagation of submerged macrophytes at the coring position (12-m water depth at present).
We interpret the enhanced terrestrial contribution (lower P aq ) at ca. 31-30, 28-26, and 17-14 kyr BP as a key to understand the occurrences of lake-level highstands during those intervals. Since vegetation in hyperarid environments reacts sensitively to slight changes in humidity, we suggest that enhanced summer humidity during those intervals caused the higher influx of terrestrial n-alkanes, which occurred synchronous to higher pollen counts and a change from alpine desert to steppe vegetation (Heinecke et al., 2018). While enhanced discharge from glacial meltwater is less likely (glaciers in the Pamirs were still either at a maximum or even expanding in the late MIS 3/early MIS 2 and in the late glacial; Dortch et al., 2013;Komatsu & Tsukamoto, 2015), an explanation of the higher lake levels are shifts in the hydrological balance toward higher precipitation and/or lower lake-water evaporation.

Changing Seasonality of the PME as Major Driver of Hydrological Changes
The simulated differences in mean annual δD-precip between the LGM and PD were approximately +11‰, according to the outputs of CAM3-iso and ECHAM5-wiso models (Table S1). This is in range of the observed offset in δD-aq (representing annual δD; Figure 2p) in the Karakul record. The +60‰ shift in δD-aq between 17 and 14 kyr BP, however, is not reflected in CAM3-iso model simulations for Heinrich event 1 (H1). Similar shifts around 30, 24, and 16 kyr BP were observed in speleothem δ 18 O records from East Asia and were hypothesized to be related to events H3 to H1 (Wang et al., 2001). Through a modeling study, Pausata et al. (2011) showed that the weakening of the Indian summer monsoon during those events leads to higher δD-precip over India and transport of anomalously D-enriched vapor to East Asia. This explanation refers to a source isotopic signature of summer precipitation and cannot be directly transferred to our study area, neither mechanistically nor with respect to the Karakul isotopic record, which lacks such a pattern in δD-terr. We hypothesize that episodes of increased summer humidity and higher lake levels at Lake Karakul were indeed a local response to northern hemispheric rapid climate events, but with more regional causes for positive excursions in aquatic isotopes.
A change of nC 23 sources (Potamogetonaceae to Characeae) could explain at least part of the large isotopic change, although limited data and contradictory results of previous studies make a full evaluation of the potential effect difficult (Liu, Yang, Cao, Leng, & Liu, 2018;. A synchronous shift toward higher values was also observed in δ 18 O of aragonite crystals (δ 18 O-ara; Heinecke, Epp, et al., 2017;Figure 2o). Between ca. 14-12 kyr BP, these two proxies for isotopic composition of lake water drift apart, visible in a sharp drop of δD-aq in contrast to only slightly decreasing δ 18 O-ara. A shrinking lake level and increase of salinity could be a reason for this pattern, as salinity is implicated to influence isotopic fractionation of aquatic macrophyte biomarkers . However, this effect is not yet fully understood, and it needs to be noted that δD-aq follows the major drop of δD-terr; hence this most likely reflects the source signal of precipitation.
We propose that changes in seasonality of precipitation provide the most likely mechanism for positive shifts in δD-aq even though it cannot be excluded that changes in aquatic macrophyte sources and salinity contributed to the signal. Isotope-enabled models and data from the Online Isotope Precipitation Calculator (Bowen & Revenaugh, 2003) reveal a strong seasonal amplitude of δD values, with summer precipitation up to 100‰ more positive than winter precipitation (Figure 3). Hence, increased summer precipitation and/or lower influx from winter precipitation best explain the +60‰ shifts in δD-aq around 31-30, 28-26, and 17-14 kyr BP (Figure 2p).
The models give further support that regional hydrology at Lake Karakul was mainly driven by changes in seasonality of effective humidity. They all indicate slightly wetter summers (increased precipitation and relative humidity) and more arid and windier winters during LGM and H1 compared to the Holocene (Figures 3  and 4b), while the annual precipitation amounts remained relatively constant (Table S1). The EC-Earth model suggests that strongly decreased evaporation during cold episodes is a relevant influencing factor behind the increased summer humidity as simulated by the models (Figure 3; Table S1). Proxy evidence points toward late glacial lake transgression starting during cold episodes, possibly as early as the LGM. The lake level highstands during the glacial period cannot be attributed to stadials or interstadials with absolute certainty due to the short-term nature of these events and uncertainties of the age-model during that period. Nevertheless, we hypothesize that similar mechanisms as in the late glacial (decreased evaporation during summers) also triggered late MIS 3/early MIS 2 lake transgressions.
On a spatial scale, most of the Central Asian mountains (i.e., the Pamirs, Hindukush, and Tien Shan) show enhanced summer precipitation during the glacial compared to the Holocene, except the eastern Pamirs where this pattern is less pronounced due to the diminished influence of the westerlies (Figures 4b and S4). In contrast, PME is clearly higher during the LGM than during the Holocene throughout whole Central Asia (Figures 3 and 4c), highlighting the relevance of this parameter as major impacting factor on transregional hydroclimate.

Conclusions
Δ terr-aq and P aq as integrative proxies for vapor source/winter westerlies intensity and catchment ecology/summer wetness, show distinct insolation-driven trends. δD-terr represents spring/summer δDprecip and its trend resembles those from other isotopic records from Asia, illustrating linkages between vapor-source changes over a large area. Based on our data it remains possible that the Indian summer monsoon touched the study area during early Holocene intensification. More likely, an intensified monsoon influenced the pathways and strengths of the westerlies, causing the most negative δD-terr values between 12 and 7 kyr BP.
δD-aq represents an integrated signal of the whole lake catchment. The large isotopic variability on glacialinterglacial scales can only be explained by significant changes in seasonality during the studied time period. Here, positive shifts at 31-30, 28-26, and 17-14 kyr BP are interpreted to be triggered by increased summer precipitation and/or low summer evaporation and decreased input from high-altitude meltwater. Synchronous higher concentrations of terrestrial biomarkers (lower P aq ) and pollen counts also give evidence for expansion of vegetation which reacts sensitively to slight increases in moisture availability in a hyperarid environment. These seasonality changes are in agreement with simulated data from CAM3-iso, ECHAM5-wiso and EC-Earth models and possibly constitute a local response to northern hemispheric rapid climate events.
Higher lake levels during these three intervals are supported by multiple indicators. For the late glacial, proxy evidence clearly indicate lake transgression during cold intervals, possibly starting as early as during LGM, long before a 14.7-kyr Bølling-warming and while glaciers were still expanding. Climate models point out decreased evaporation as a major factor that increases summer effective moisture availability and profoundly influences the hydrological states of lakes.

Data Availability Statement
Data presented in this study are available on PANGAEA: https://doi.pangaea.de/10.1594/PANGAEA.907737 from the Swedish Research Council VR project 2013-06476. The simulations with CAM3 and EC-Earth were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at NSC and Cray XC30 HPC systems at ECMWF. We acknowledge two anonymous reviewers whose comments improved the manuscript.

Erratum
In the originally published version of this article, coauthor Ilhomjon Rajabov's name was spelled incorrectly. This error has since been corrected, and the present version may be considered the authoritative version of record.