What are the drivers of Caspian Sea level variation during the late Quaternary?

Quaternary Caspian Sea level variations depended on geophysical processes (affecting the opening and closing of gateways and basin size/shape) and hydro-climatological processes (affecting water balance)


Introduction
Understanding drivers and future trajectories of sea and lake level change is of paramount significance for mitigating future socio-economic disasters and for planning conservation measures. Long-term changes are partly related to geophysical processes altering the shape and size of sea/lake basins. Changes to the volume of water held within those basins also impact sea/lake level. Volume changes result from either density-driven variation (salinity and/or temperature changes), or from changes in mass (e.g., due to variation in water locked up in ice-sheets and glaciers). Human interference with water storage on land for agriculture, industry, and direct human consumption also has the potential to contribute to current sea level change (Pokhrel et al., 2012).
The Caspian Sea is an endorheic long-lived lake with sea-like features (saline, marine like faunas), whose magnitude of water level variations has been larger and of higher frequency compared to the marine realm (Kroonenberg et al., 1997(Kroonenberg et al., , 2000 (Fig. 1). Large populations are reliant on the water resource within the drainage basin for irrigation, industry, and domestic uses (Koriche et al., 2021a). Understanding the natural and anthropogenic controls on Caspian Sea level is therefore of critical importance. Using past variations in sea level can help understand drivers (processes and settings) that affect Caspian Sea level dynamics. However, scarcity/ ambiguity of age and sea level data in the Pontocaspian region (Black Sea and Caspian Sea basins) have made it difficult to identify consistent periods of sea level change in the palaeo-environmental record (Krijgsman et al., 2019).
During the Quaternary, the record suggests there were typically extended periods of isolation and short intervals of connection between lake basins in the Pontocaspian domain (Bezrodnykh et al., 2004;Rychagov, 1997a;Svitoch, 2013). In particular, the Caspian Sea experienced large variations in water level from tens to hundreds of meters on various time scales (Bezrodnykh et al., 2004;Fig. 1. Adapted from Fig. 6 of Krijgsman et al. (2019). a) Black Sea and b) Caspian Sea level reconstruction from Pleistocene to Holocene. Water flux between Caspian Sea and Black Sea are indicated by blue arrows facing to the left (or both ways for bi-directional flow). The Caspian Sea level curves are reconstructed based on various palaeo-environmental records. The uncertainty may depend on the dating method employed and the bias in the quantification of the Caspian Sea level variation. For example, a broad array of fauna and flora are used as sea-level indicators from coastal regions (e.g., corals, benthic foraminifers, ostracods, diatoms and mangroves) and from marine terraces (e.g., sedimentary and or stratigraphic features), which are dependent on the nature of the indicators, methods used, palaeo period considered and the regional environmental setup. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) Svitoch, 2010;Varuschenko et al., 1987;Yanina, 2014) (Fig. 1). Several prolonged Caspian Sea highstand phases (e.g. Apsheronian, Bakunian, Khazarian, and Khvalynian) have been distinguished in the Quaternary (Forte and Cowgill, 2013;Krijgsman et al., 2019;Kroonenberg et al., 1997;Yanina, 2012).
Caspian Sea level variation is unique in nature and its evolution is complex, due to the various driving processes and the environmental setting as part of a closed basin. Understanding the relative contributions of the different drivers of Caspian Sea level change requires a multi-disciplinary approach (e.g., hydroclimate modelling as well as geological and geomorphological analysis). A considerable number of studies have explored various aspects of Caspian Sea changes (e.g., Arpe et al., 2014;Arslanov et al., 2016;Kislov and Toropov, 2007;Kislov et al., 2014;Rychagov, 1997a;Tudryn et al., 2013;Yanina, 2014). However, there are conflicting ideas concerning the relative importance of different factors driving past Caspian Sea level variability (discussed in Section 2), as well as the mechanisms by which they impacted basin connectivity and Caspian Sea level variation. Below, we outline plausible sources of runoff and flow directions (see Fig. 2; sources are indicated by the numbers in a square brackets) that potentially played a role in controlling Caspian Sea level variation during the Quaternary (especially during glacial/deglacial periods): meltwater from the Eurasian (Fennoscandian) Ice-sheet (Fig. 2, [1]); river basin and drainage reorganization due to ice-sheet damming and accumulation; expansion of river basin area due to topographic change produced as a result of ice-sheet loading/unloading; formation of closed depressions and proglacial lakes outburst (Fig. 2, [2] and [4]); glacial melting (Fig. 2,[3]) and hydroclimate change in the neighbouring Himalayan and Amu-Darya basins, that were likely hydrologically connected to the Caspian basin at various times; sedimentation and structural change of spillways due to isostatic rebound (e.g. Manych-Kerch) (Fig. 2 [5] and [6]); Any of the above changes would have been in addition to the hydroclimate variation (precipitation-evaporation balance) controlled by glacial-interglacial cycles and millennial scale forcing. We anticipate that Caspian Sea level variations originate from the sources of runoff and changes in drainage extent outlined in the list above. To disentangle these factors and understand their relative importance, we firstly performed a synthesis of palaeo lake level data covering the last 25 kyr (see section 4.1). Secondly, we examined the impacts of hydro-climatological changes using palaeoclimate model simulations (see section 4.2). Thirdly, we evaluated the impact of the Fennoscandian ice-sheet on catchment dynamics and meltwater contributions (see section 4.3). Finally, we constructed a model of Caspian Sea level in order to quantify the cumulative impact of the different drivers (see section 4.4).

Previous studies
As previously indicated, Caspian Sea level has varied during the Quaternary due to both geophysical processes (that have opened and closed gateways) and hydroclimate change induced by glacialinterglacial cycles (Badertscher et al., 2011;Kroonenberg et al., 2005;Richards et al., 2017;Rodionov, 1994;Yanina et al., 2018;Yanina, 2014). The gateway that controls the connection between the Ponto-Caspian seas is the Manych-Kerch (Fig. 2). This gateway is affected by tectonics (Le Pichon et al., 2016;Svitoch and Makshaev, 2011), deposition, and erosion (e.g., Svitoch and Makshaev, 2011); these mechanisms modify topography and therefore may lead to changes in drainage basin size or shape, or to changes in lake bathymetry. The architecture of the gateways impacts the base-level of upstream basins directly. The present rates of the tectonic movements are small (<1 mm per year) over most parts of the Manych-Kerch straits/depressions (Svitoch, 2013). However, in earlier time periods (e.g. the early Pliocene), geophysical processes (sedimentation and subsidence) played a vital role in Caspian Sea level variation (Kroonenberg et al., 2005;Nadirov et al., 1997;Richards et al., 2017;Vincent et al., 2010).
Many studies suggest that, hydroclimatic processes play a significant role in controlling Caspian Sea level variation (e.g., Badertscher et al., 2011;Dolukhanov et al., 2010;Jorissen et al., 2020;Kislov and Toropov, 2007;Kislov et al., 2014;Yanina, 2012;Yanko-Hombach and Kislov, 2018;Mayev, 2010;Levchenko and Roslyakov, 2010). The regional water budget is mainly influenced by freshwater balance, the relationship between evaporation and lake area, melting of glaciers and/or ice-sheets, and permafrost melting (Kondratjeva et al., 1993;Kroonenberg, 2012;Rodionov, 1994;Yanina, 2012). Present-day Caspian Sea level variation is mostly controlled by contributions from more than a hundred rivers that discharge into the Caspian, although close to 80% of the discharge originates from the Volga (Arpe et al., 2012;Rodionov, 1994). The discharge from these rivers is affected by hydroclimate processes at global and regional scales. Increased river discharge has been reported for interglacial periods, as well as during melting of ice sheets in periods of deglaciation (Krijgsman et al., 2019 and references therein). Other additional sources of runoff during the Quaternary include previously hydrologically connected neighbouring river basins (e.g., Amu-Syr-Darya basin; Fig. 2) that are not connected in the present day (Kroonenberg et al., 2005;Leroy et al., 2013;Torres et al., 2007). For example, runoff to the Caspian Sea via the Uzboy River from the Aral Sea and Amu-Darya occurred briefly at around 5.8 to 5.2 and 4.5 to 3.8 14 C kyr BP (Leroy et al., 2013 and references therein).
According to Mangerud et al. (2004), the Caspian Sea received significant discharge from ice-dammed proglacial lake outbursts from the south of the Siberian (BarentseKara) ice-sheet via the Aral Sea during 90e80 kyr BP and around 18e17 kyr BP. Svitoch (2009) alternatively suggested that the Khvalynian transgression was never as a result of ice-dammed proglacial lake outbursts, but fed by discharge from Amu-Syr Darya drainage basin through Lake Sarykamysh (located approximately midway between the Caspian Sea and the Aral Sea). He also points that, this runoff source was not significant and had only a small contribution to the Khvalynian transgression. This is also supported by the review by Panin et al. (2020) on the routes of runoff to Caspian Sea. They also suggest that the amount of runoff that joined Caspian Sea via Lake Sarykamysh (Uzboi valley) would have only raised Caspian to À22 m above mean sea level, which corresponds to the maximum Holocene Caspian Sea level, but cannot explain the Pleistocene Caspian Sea level changes, indicating that this route can be ruled out as the source of runoff for the Khvalynian transgressions. Tudryn et al. (2016) indicated that the Caspian Sea acted as a final trap for south-eastern Scandinavian ice-sheet meltwaters during the last deglaciation. The meltwater discharge joined the Caspian Sea via the Volga River from the LGM until~13.8 cal kyr BP, as suggested by Tudryn et al. (2016) and Soulet et al. (2013). Most of the late Pleistocene and Holocene water level fluctuations of Caspian Sea are thought to be due to hydroclimate change (Krijgsman et al., 2019;Yanchilina et al., 2019;Yanina et al., 2018). This is also supported by various palaeoclimate model studies (Kislov and Toropov, 2007;Renssen et al., 2007). Runoff from the catchment was found to be the primary driver of Caspian Sea level variation during the mid-Holocene and LGM, based on Palaeoclimate Model Intercomparison Project 3 (PMIP3) models (Kislov and Toropov, 2007).
In summary, although a considerable number of studies have focussed on exploring the main source of Caspian Sea level variation, it is likely that there have been a number of contributing mechanisms. In this study, we build on previous work to consider the cumulative impacts of the various potential drivers of Caspian Sea level variation and changes in their relative contributions through the last deglaciation.

Study area
The Caspian Sea formed as an isolated closed basin in the early Pliocene (Popov et al., 2006;van Baak et al., 2016). It is the largest inland water body in the world and has a present-day catchment area of approximately 3.6 Mkm 2 (see Fig. 2 key palaeo and presentday features). The vast catchment covers a number of climatic zones, from a temperate continental climatic zone over the Volga basin, to moderately warm climate over the western area, subtropical humid climate over the southwestern and the southern regions, and semi-arid climate over the eastern part (Chen and Chen, 2013;Kosarev, 2005;Leroy et al., 2007Leroy et al., , 2020). Currently, Caspian Sea level is~28 m below mean sea level. The Caspian system is very sensitive to changes in evaporation over the Caspian Sea (Arpe et al., 2012;Chen et al., 2017;Rodionov, 1994). Due to the shape of the lake bathymetry, as Caspian Sea level changes (between À100 m and þ50 m above mean sea level), its surface area changes significantly and non-linearly (Fig. 3). This has important implications for determining the amount of evaporation from the Caspian Sea (a larger surface area often corresponds to strongly increased evaporation), as well as climate feedbacks with the atmosphere. Climate impacts occur in both the local-scale atmospheric water-balance and large-scale atmospheric circulation patterns (Koriche et al., 2021b). The relationship between Caspian Sea volume, surface area, and water level is used within the hydrological model developed and outlined in the next section.

Research methods
In this section we present the datasets used, hydrological modelling scheme implemented, and methods used for data analysis. The first part of the study focuses on a synthesis of Caspian Sea level data from the literature. The second part involves analysis of hydroclimate evolution over the Caspian Sea basin based on climate model simulations of the last glacial cycle (120e0 kyr BP). The third section looks at the impact of the Fennoscandian ice-sheet on the re-organization of the Caspian Sea drainage system and its contribution of meltwater to the basin over the last 25 kyr. The fourth part describes a hydrological model that we constructed to explore the cumulative impact of hydroclimate change, meltwater input, and river drainage system rearrangement on Caspian Sea level. The potential contribution due to ice-sheet meltwater via Aral Sea and the east of the Caspian Sea were not examined in this study as the topographic information is only based on coarse (1 ) resolution, which limits the routing of runoff over the palaeo period.

Reconstruction of palaeo-Caspian Sea level
We reconstructed Caspian Sea level over the last glacial cycle based on the limited available palaeo-environmental records from published studies. Table 1 contains details of the data sources and methods used to infer past Caspian Sea level as well as stratigraphic time data. We used the original chronologies associated with each dataset. Data points with age uncertainties of ± 5 kyr or larger were excluded, though these data sources are included in the Supplementary information, Table S1, for completeness. Consequently, only Caspian Sea level reconstruction points covering the last 25 kyr were suitable for inclusion in the final compiled dataset, presented in section 4.1.

Palaeoclimate modelling
Climate model data was extracted from snapshot simulations covering the last 120 kyr BP (Davies- Barnard et al., 2017;Singarayer and Valdes, 2010;Singarayer et al., 2017) performed using the Hadley Centre coupled atmosphere-ocean-vegetation climate model, HadCM3 (Gordon et al., 2000;Pope et al., 2000). The simulations included the MOSES2.1 (Met Office Surface Exchange Scheme) land surface scheme and the TRIFFID (Top-down Representation of Interactive Foliage and Flora Including Dynamics) dynamic vegetation model (Cox et al., 1999). HadCM3 has 19 vertical levels with a horizontal resolution of 3.75 Â 2.5 in the atmosphere and 20 vertical levels with a resolution of 1.25 Â 1.25 in the ocean. The model experiment produced 65 climate snapshots over the last 120 kyr BP, as producing a global time series using a fully coupled complex climate model transient simulation is far too computationally expensive. Each was forced with boundary conditions pertinent to the time period being simulated, including insolation seasonality, atmospheric greenhouse gas concentrations based on ice-core data, and ice sheets. From 120 to 80 kyr BP, simulations were performed at intervals of 4 kyr; between 80 and 22 kyr BP simulations were performed at 2 kyr intervals; and then simulations were created for every 1000 years from 22 to 0 kyr BP (Singarayer and Valdes, 2010;Singarayer et al., 2017). The key assumption in snapshot experiments is that the climate is in equilibrium with the boundary conditions, and that the final climate is largely independent of the initial conditions. Both of these assumptions can be challenged during the late Pleistocene period (last deglaciation period) as the time is characterized by rapid climate change like the Heinrich event and Younger Dryas where the changes were happening in a few years (Steffensen et al., 2008). However, the use of freshwater hosing simulations for these particular periods (where equilibrium is not reached before taking the mean climatologies) somewhat corrects for this. In addition, we incorporated millennial-scale forcing of freshwater into the North Atlantic to simulate iceberg influx during Heinrich event 2, Heinrich event 1, and the Younger Dryas. During such millennial-scale Heinrich events, palaeodata indicate that the Atlantic Meridional Overturning Circulation (AMOC) collapsed (Seidov and Maslin, 1999), resulting in cooler temperatures and drier conditions over Europe and much of the northern hemisphere (e.g., S anchez-Goñi et al., 2000) over an extended period of time (>1000 years). Such events therefore had the potential to influence the surface water balance over the Caspian basin due to changes in precipitation and evaporation, and as a result, impact Caspian Sea level. In HadCM3 the influx of iceberg discharge is simulated as a negative salinity flux equivalent to the addition of 0.3 Sv over the surface of the North Atlantic (between 50 and 70 N). The result of this is a reduction in AMOC overturning strength. A limitation of this technique is that it assumes that all millennial scale variability is driven by changes in the AMOC. The decrease in AMOC strength varies depending on the background climate state. For the Younger Dryas simulation, the AMOC experiences a reduction from 18 Sv to 5 Sv on average.
The glacial cycle HadCM3 simulations have been evaluated by model-data comparisons in a number of studies covering various climate fields and geographical regions (e.g., Hoogakker et al., 2016;Singarayer and Burrough, 2015;Singarayer and Valdes, 2010;Singarayer et al., 2017). The modelled global average glacialinterglacial temperature change in this version of HadCM3 is 5.5 C, which is in the middle of the range inferred from palaeodata (Masson-Delmotte et al., 2013). Polar temperature trends compare well to ice core records, although the magnitude of change is underestimated (Singarayer and Valdes, 2010), as it is in most climate models. Temporal variations in large-scale circulation such as the position of the intertropical convergence zone and monsoon intensity at locations from East Africa to China are well simulated by HadCM3 (Singarayer and Burrough, 2015;Singarayer et al., 2017). The glacial-interglacial changes in climate have been used to drive biome simulations that match well with biome reconstructions from pollen over Europe and Asia (Hoogakker et al., 2016). Comparison with modern observations is described in . HadCM3 compares well with newer state-ofthe-art models and is many times faster to run, therefore making it useful for simulating palaeoclimate scenarios over multiple time periods.
Note that HadCM3 does not include interactive ice sheets or carbon cycle. Thus, the ice-sheet evolution and greenhouse gases have to be prescribed (Singarayer and Valdes, 2010). In these simulations (with HadCM3), the ice sheet is based on ICE-5G (Peltier, 2004). Prescribing ice-sheet evolution helps to appropriately represent ice-sheet extent and palaeo-topography. However, representing the actual ice-sheet meltwater and its contribution to the basin water budget is vital and is not included in the HadCM3 climate simulations. Therefore, we separately estimated the meltwater contribution to Caspian Sea level using ice-sheet reconstructions from ICE-6G_C (Argus et al., 2014;Peltier et al., 2015).
We carried out a model-data comparison of the surface water balance in the present day to assess whether it was appropriate to incorporate a bias correction into our analysis. Prior to bias correction, the model P-E (precipitation minus evaporation) outputs were downscaled from the standard HadCM3 resolution (3.75 Â 2.5 ) to 6-arcminutes by a first-order conservative interpolation method (remapcon) using the Climate-Data-Operator (CDO) software. The conservative interpolation method works well for flux conservation and interpolation to higher resolution  Kroonenberg et al. (2008) [g] Multiple sources* Early Holocene: 12 to 7.5 (Mangyshlak regression) (Bezrodnykh et al., 2004(Bezrodnykh et al., , 2020Bezrodnykh and Sorokin, 2016;Kroonenberg et al., 2005Kroonenberg et al., , 2008Sorokin, 2011;Tudryn et al., 2013) [h] Multi-proxy (foraminifera, Diatoms, organic matter) Coring, sedimentological analysis and bio-facies 14 C LGM and mid-Holocene Kakroodi et al. (2015) [i] Marine shells Analysis of the geology and geomorphology terraces; levelling of terraces and related shore-lines, 14 C Holocene Rychagov (1997a) [  (Jones, 1999). Three different sources of data were then used to estimate observational surface water balance over the Caspian basin. These include: (1) the 20th century observed mean Caspian Sea level change record (where average year-to-year variations in Caspian Sea level give the annual water balance) combined with average human water use from during the 20th century based on the work of (Rodionov, 1994) and (Shiklomanov, 1981); (2) the ECMWF 20th century reanalysis (Poli et al., 2016); and (3) the 5th generation ECMWF reanalysis data (ERA5; Copernicus Climate Change Service (C3S), 2017; Hersbach et al., 2020). By using these three observational datasets for bias correction we were able to quantify a measure of uncertainty in the observed annual water budget for the Caspian basin. Results are described in Section 4.2.
Since there are no direct observational datasets to use for bias correction of the palaeo-time slices, we computed the mean annual water balance (P-E) anomaly by subtracting the pre-industrial P-E from the P-E at each palaeo-time slice. The assumption is that, the difference between the present-day mean P-E (~2950 m 3 s -1 ) and pre-industrial HadCM3 P-E (~3300 m 3 s -1 ) is small. Therefore, the observed present-day water balance was added to each of the palaeo P-E anomalies. As indicated above, due to shortage of observed palaeo datasets, this approach assumes that the bias is constant throughout all palaeo-time slices, however, this is not necessarily the case in the reality, as the climate over this period varied due to differences in the forcing and boundary conditions. We compare results with and without bias correction. The results are discussed further in Section 4.4.

Palaeo drainage basin characterization
Ice sheets play an important role in reshaping the river drainage system. They also contribute vast quantities of meltwater during deglaciation periods, resulting in abrupt sea level and climate change. We evaluated the impact of the Fennoscandian ice-sheet on the dynamics of Caspian Sea drainage area evolution and its meltwater contribution. The drainage basin and the associated palaeo-rivers were reconstructed by combining a high resolution (0.5-arcminutes) digital elevation model (DEM) from the General Bathymetric Chart of the Oceans (GEBCO) gridded merged topographic and bathymetric data set (GEBCO, 2019) and palaeo icesheet reconstruction based on postglacial rebound model ICE-6G_C (Argus et al., 2014;Peltier et al., 2015). ICE-6G_C is a global ice-sheet and glacial isostatic adjustment (GIA) reconstruction created by combining geological data and geophysical model (Argus et al., 2014;Peltier et al., 2015). It has shown better performance compared to its predecessors (ICE-3G, ICE-5G) and closer to glacial geological data over the Laurentide icesheet (Wickert, 2016). Over the Fennoscandian, the ICE-6G_C total volume is 15% less than ICE-5G during 21 kyr BP and early deglaciation period ( Fig. 4a). As shown in Fig. 4b, the magnitude of meltwater contributions differ during the LGM and early deglaciation period. Though the difference is smaller, it could imply that the simulated runoff by HadCM3 would be different for both ice-sheet reconstructions. However, in HadCM3, the ice-sheet extent was prescribed, and interactive ice sheet evolutions were not considered. Therefore, the impact on the simulated runoff would be less significant.
Two sets of palaeo-topography were created based on 1) including ice-sheet thickness without GIA, and 2) including the topographic anomalies from present day where both ICE-6G_C icesheet thickness and GIA were considered. Since the reconstructed ice-sheet geometry was only available at a relatively coarse resolution (1 x 1 ), combination with high-resolution elevation data was necessary as some topographical features cannot otherwise be sufficiently well represented, such as the valleys of major river systems. High-resolution elevation data is only available from present-day observations. In order to combine both data sets, the coarse resolution ice-sheet palaeo time slices were interpolated using an iterative (7 steps) bi-linear method to help remove stepwise discontinuities that may introduce unrealistic flow direction during drainage basin creation (Wickert, 2016). High-resolution anomalies of palaeo-time slices minus pre-industrial time slices were then added on to the high-resolution modern topographic data to create the 0.5 arc-minute resolution palaeo-topographic maps (Fig. 5).
The characterization of the drainage of the Caspian Sea basin was performed by employing Geographic Resources Analysis Support System (GRASS) software (Metz et al., 2011;Neteler et al., 2012) and using the high resolution palaeo-topographies as key inputs. GRASS is an open-source geospatial analysis software recommended for characterizing hydrologically linked drainage basins. From this analysis, drainage basin area and the associated river networks covering the period from 25 kyr to present were created. The derived drainage extent was used to perform analysis of palaeo-water budget and ice-sheet meltwater contributions.

Caspian Sea level modelling
Caspian Sea volume and level were simulated by integrating all the plausible sources of runoff due to palaeoclimate change, drainage river system rearrangement, and ice-sheet meltwater contributions to the water budget. Here, the Caspian Sea model (Eqs. (1) and (2)) is based on fluxes of precipitation minus evaporation (P-E) over the Caspian Sea drainage basin and meltwater from the ice sheet within the drainage basin. However, abrupt climate change and/or high discharges of meltwater can increase Caspian Sea level enough that water starts to spill over to the hydrologically linked neighbouring basin (the Black Sea) via the Manych-Kerch spillway. Equally, there is a lower volume limit set of 0 m 3 , which may be reached over time given sufficiently negative fluxes. These limits are expressed in Equation (2).
where: CSV is Caspian Sea volume, P is precipitation, E is evaporation (from canopy, soil and water body), A is surface area over the land or sea part of the basin, M is the ice-sheet meltwater flux, B is a bias correction based on modern day observations, Dt is the time step, and CS_max is the maximum volume in the Caspian Sea before overspill occurs. Once CSV is calculated, the Caspian Sea level at each time step is calculated by interpolation of the volume-sea level curve in Fig. 3. As mentioned previously, the area of the Caspian Sea varies considerably with changing Caspian Sea level (Fig. 3). Any changes in sea surface area would significantly impact the total evaporation per time step. Accordingly, we use a variable area (calculated from the volume at the current time step) and multiply the water balance (P-E) over the Caspian Sea by the surface area (at the previous time step; Eq. (1)), and similar for the land area water balance. In order to avoid numerical instabilities, we used a time step (Dt) of 20 years. As the simulated climate and ice-sheet states are averages at 1000-year intervals, the precipitation, evaporation, and ice-sheet meltwater fluxes were linearly interpolated to 20-year intervals.
Previous studies have shown that the contribution of groundwater to the Caspian Sea in the late Quaternary is estimated to be small (Golovanova, 2015;Zekster, 1995). Therefore, we assume that there is no significant groundwater reservoir on these timescales and have not included a groundwater component in the water balance analysis.

Synthesis of palaeo lake level data
As previously outlined, present day Caspian Sea level is at an elevation of~28 m below mean sea level and the water balance is controlled mainly by inflow from major rivers, such as the Volga, Kura, and Ural, balanced by evaporation over the Caspian Sea (Arpe et al., 2014;Chen et al., 2017;Rodionov, 1994). Over the last thousand years the Caspian Sea level varied between 19 and 28 m below mean sea level, based on Caspian Sea level reconstruction from combined historical documents and geological records by (Naderi Beni et al., 2013) (Table 1). The rate of Caspian Sea level change earlier in the Holocene was much higher compared to present-day Caspian Sea level change (variation of around 3 m in the 20th Century). The deepest lowstand recorded (referred to as the Mangyshlak regression) for the Caspian Sea during our time frame of interest was during the early Holocene (Fig. 6). This lowstand was approximately 80e113 m below mean sea level (Kroonenberg et al., 2008). The timing of the Mangyshlak lowstand varies between studies from 12 kyr BP to 8 kyr BP (e.g. Bezrodnykh and Sorokin, 2016;Yanina, 2014;Kroonenberg, 2012) (Fig. 6). For instance, Bezrodynkh and Sorokin (2016) indicated that the start of the Mangyshlak regression was at around 12.4 kyr with a level of around À40 m which then dropped rapidly by 40e44 m over 700e1000 years, and they indicate an end of the regression around 9.5 kyr BP (see Table S1 for other sources). However, Rychagov (1997b) suggests that maximum lowstand was around 12 kyr BP. The second lowest Caspian Sea level recorded in the Holocene was during the warm medieval period (at À42 m a.s.l), around 1.4 kyr BP (Kakroodi et al., 2012). Numerous Caspian Sea highstands of up to 25 m below mean sea level were also recorded semi-regularly throughout the Holocene (Kakroodi et al., 2012;Rychagov, 1997a).
A series of palaeo-geographic events took place during late Pleistocene, including the Late Khazarian, early Khvalynian and late Khvalynian transgressive stages, and the Atelian, Enotaevka and Mangyshlak regression stages (see Fig. 1 and Table S1 for suggested age and the corresponding Caspian Sea level). During these periods, the Caspian Sea experienced extreme water level change, ranging from þ50 m to~-130 m above mean sea level. Various studies explored timing and magnitude of Caspian Sea level events from depositional records (e.g., Bezrodnykh et al., 2015;Bezrodnykh et al., 2004;Mamedov, 1997;Shkatova, 2010;Sorokin, 2011;Tudryn et al., 2013;Yanina et al., 2018;Yanina, 2014) but no consensus exists over the Caspian Sea level curves. The ages of the Atelian regression and the early Khvalynian transgression are still controversial (Krijgsman et al., 2019 and references therein;Yanina et al., 2018). Many authors (e.g., Dolukhanov et al., 2010;Kroonenberg et al., 1997;Tudryn et al., 2013;Yanina, 2014) reported the existence of a þ50 m highstand of the Caspian Sea for the early Khvalynian transgression with age estimates ranging between >40 kyr BP to 17 kyr BP (see Table S1).
In summary, the Caspian Sea level variation during the late Pleistocene glacial period was much larger than the Holocene, with the potential of overspill periods to the Black Sea basin (Fig. 6). The variation reached more than 70 m of water level change. The extreme (Mangyshlak) lowstand at the Pleistocene-Holocene transition suggests that the climate rapidly became drier. This was then followed by transgressive and regressive phases ranging from À19 to À41 m a.s.l., controlled by hydroclimate change over the Caspian Sea basin (Bezrodnykh et al., 2020;Kroonenberg et al., 2008) combined with intermittent contributions from hydrologically linked basins such as the Amu-Syr Darya river basin around 4 14 C kyr BP .

Hydroclimate analysis based on HadCM3 simulations
Until recently, palaeoclimate model simulations were typically focused on snapshot simulations of a few well-studied periods, such as the LGM and Mid-Holocene, which have been the focus of international modelling efforts as part of the PMIP collaboration (Braconnot et al., 2012). However, here we want to understand the cumulative impact of changing climate conditions over the late Quaternary. In order to advance our understanding of the impact of climate change on Caspian Sea level evolution, we have analysed the basin water budget variation during the last 25 kyr using simulations from a global coupled ocean-atmosphere-vegetation climate model (HadCM3), based on 25 time-slice simulations (as described in section 3). The modelled annual mean water budget variations presented in Fig. 7a (black sold line) are integrated over the present-day Caspian Sea basin.
We observe that the pre-industrial (0 kyr BP) model estimate of P-E over the present Caspian Sea basin is higher than the rest of the time slices considered. The P-E balance is generally positive in interglacial periods (see Fig. S1 for model simulation results covering the last 120 kyr) and then lower during the glacial period (Fig. 7a). The spatial pattern (Fig. S2) of precipitation and evaporation change is somewhat heterogeneous and varies between time slices, but there are common features. In general, both precipitation and evaporation over land in the drainage basin decrease during glacial times compared with the pre-industrial. However, precipitation declines to a greater extent than evaporation, resulting in drier conditions (an overall decrease in P-E). Over the northern Caspian Sea itself and surrounding land, both precipitation and evaporation increase in the modelled glacial time slices compared to the pre-industrial, primarily as result of an expansion of the prescribed lake grid cells. However, in this region the increase in evaporation is larger than that of precipitation due to the inundation of the land by the Caspian Sea, and so here also the change in P-E is negative (i.e. drier conditions). Over the southern lake grid cells glacial P-E is higher than pre-industrial. Here, there is no expansion of the lake to replace land grid cells, and while precipitation decreases, the evaporation decrease is greater due to reduced glacialstage temperatures (Fig. 7a). The rapid climate change events that are simulated (Heinrich events 1 and 2, and the Younger Dryas) produce drier conditions (negative change in P-E) that is roughly half of the total glacialinterglacial amplitude, but on a shorter timescale. The freshwater  Table 1 and listed below). The error bars indicate the uncertainty in the reported age records except for 'Enotaevkan' regressions, where various studies reported different lowstand timings, which are combined into a single uncertainty bar for age and Caspian Sea level. The grey-shaded bar indicates the period 'Mangyshlak' regression lasted based on various studies (see Table S1). The pink-shaded bar indicates Caspian Sea overflow into the Black Sea basin according to Badertscher et al. (2011)  flux forcing introduced in the North Atlantic induces a rapid collapse of the AMOC and reduced meridional heat transport, which produces widespread temperature decreases throughout the northern hemisphere. The temperature reductions over Europe and the Caspian basin result in decreases in both evaporation and precipitation. As with the multi-millennial scale changes, the reduction in evaporation is not as large as precipitation, producing an overall drying (i.e. negative change in P-E). Despite each of the simulated millennial events being forced with the same magnitude of freshwater flux, the magnitude of the P-E impacts varies. This is likely due to a dependence of the response on the background climate state (e.g., Gong et al., 2013), although further exploration of this is beyond the scope of this study.
In Fig. 7b, we present the annual mean HadCM3 P-E, adjusted using bias correction (as described in Section 3.2.4), and a plausible range of uncertainty of P-E for the present-day Caspian Sea basin. The mean annual P-E with bias correction from ERA5 (Hersbach et al., 2020) is given as the middle line in Fig. 7b (black solid  line). The upper bound of the uncertainty region (shown in grey shading) is given by the annual mean water budget (~6835 m 3 s -1 ) derived from historical Caspian Sea level records combined with an average annual rate of human water extraction (as described in Section 3.2.4). The lower bound is a plausible minimum value (~778 m 3 s -1 ) that is obtained by considering present day observed mean annual Caspian Sea level change without any human extraction (i.e. assuming that either human extraction is negligible or that all human water extraction flows directly back into the Caspian Sea). The raw pre-industrial modelled P-E is closest to the surface water budget calculated from ERA5 and sits in the middle of the uncertainty bounds. The simulated water budget based on the climate model P-E does not match with the palaeodata based Caspian Sea level graph in Fig. 6. Caspian Sea level data suggest a strong positive water budget with a period of Caspian Sea over-spill to the Black Sea during the deglaciation period (20-14 kyr BP), whereas the modelled P-E suggests that the glacial and deglacial are considerably drier than the Holocene. Of course, the Caspian Sea level record is integrating year-to-year changes in budget through time and so we would not expect the temporal evolution of the modelled P-E to produce a direct correspondence with it. However, there is no suggestion that the climate change driven changes in the water budget within the current extent of the Caspian Sea basin resulted in over-spill of the Caspian Sea in the early deglaciation phase. One interpretation of the model-data difference is that there are missing components in the modelled water budget. Potential contributions could come from (1) ice-sheet meltwater (that is neither prescribed nor modelled interactively in the climate model) and/or (2) the impact of topography changes (due to isostatic adjustment and ice-sheet configuration) on drainage system reorganization and basin area change. These issues are addressed in section 4.3.

Evaluation of Fennoscandian ice-sheet impact on the Caspian Sea
Previous studies (e.g., Wickert, 2016) have demonstrated that the North American ice-sheet complex played a vital role in restructuring rivers and drainage basin systems by altering the topography of the Earth's surface. Below we use similar methods to explore the impact of the Fennoscandian ice-sheet on reorganizing the palaeo-drainage of the Caspian Sea basin during the last deglaciation period. We explore the potential impact of the direct contribution of water as the ice melts. Detailed explanations are presented in the following sections (4.3.1 and 4.3.2).

Reconstruction of palaeo Caspian Sea basin and river flow direction
Reconstruction of the palaeo-Caspian Sea basin and river flow directions was performed using the method described in Section 3 at 1000 yr intervals over 25 kyr BP to 0 kyr BP. The spatial maps presented in Fig. 8 show the evolution of the Caspian Sea Basin and river network for key periods, including the LGM. Readers are advised to refer to the supplementary movie showing the full time series reconstruction from 25 kyr BP to present day. There was massive expansion of the drainage basin and river systems during the last glacial maximum and deglaciation (Fig. 8).
Supplementary video related to this article can be found at https://doi.org/10. 1016/j.quascirev.2022.107457 We considered two scenarios of palaeo-topography, both of which include ice-sheet topographic change but only one in which glacial isostatic adjustment is also incorporated. The spatial maps for these scenarios, presented in Fig. 8a and b, were constructed using the ice-sheet history from ICE-6G_C (Argus et al., 2014;Peltier et al., 2015;Stuhne and Peltier, 2015) combined with a present-day high-resolution digital elevation model (DEM). Fig. 8a shows spatial map of the palaeo Caspian Sea drainage system derived from topography created by combining present-day GEBCO digital elevation model (DEM) and ice-sheet thickness history from ICE-6G_C, hereafter referred as CSB-noGIA. In this scenario the effect of glacial isostatic adjustment is not considered. The spatial maps shown on Fig. 7b are derived from ICE-6G_C topography corrected for glacial isostatic adjustment (Argus et al., 2014;Peltier et al., 2015) combined with the present-day GEBCO DEM, hereafter referred as CSB-GIA. Another set of spatial maps (Fig. S3) are included in the supplementary information; these were made from the paleotopogrpahy calculated by Wickert (2016) in the course of his work to reconstruct past North American river systems. These were also prepared using ICE-6G_C (VM5a). GIA was calculated following Kendall et al. (2005), and interpolation and drainage calculations were completed following Wickert (2016). The associated drainage basins and river discharges have been published previously, as applied to paleoclimate GCM studies, by Ivanovic et al. (2017Ivanovic et al. ( , 2018aIvanovic et al. ( , 2018b, though using global-scale reconstructions that also incorporate potential inflow from highmountain Asia. This third combination produced drainage basins similar to CSB-noGIA and so is not discussed in the main text. This difference highlights the sensitivity that drainage networks in lowrelief regions may have to subtle shifts in the geometry of reconstructed ice sheets and GIA (cf. Wickert et al., 2013Wickert et al., , 2019.  (Argus et al., 2014;Peltier et al., 2015;Stuhne and Peltier, 2015); here, glacial isostatic adjustment is not considered. The second row (b) adds glacial isostatic adjustment from Peltier et al. (2015) to the topography used for our flow-routing calculations. Drainage basins are delimited in dark grey. Shaded areas and blue lines represent the respective catchment and river flow directions. The white and blue colour backgrounds depict the present-day land and Caspian Sea masks respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) In both scenarios (CSB-noGIA and CSB-GIA) the Caspian Sea basin expanded farther north of the present-day basin land area. The total area of the basin increased by 60e70% of current basin size (Fig. 9). The basin decreases in size earlier in the deglaciation when the glacial isostatic adjustment component is included. The reorganization of the Caspian Basin has implications for Caspian Sea level variation, as the drainage system experiences additional river flow and routing of meltwater into the basin with a larger drainage area.

Impact on water budget
In this section we present an evaluation of ice-sheet and topographic impacts on the palaeo-water budget of the Caspian Sea basin. Fig. 10a, shows the total volume of accumulated ice-sheet within the palaeo-Caspian Sea basin area (CSB-noGIA and CSB-GIA) and from the entire Fennoscandian ice-sheet. The substantial amount of ice-sheet accumulation during the glacial period resulted in an increase of elevation that caused the Caspian Sea area to expand (Fig. 8). North-flowing rivers were redirected southward and this led to an increase in runoff due to ice-melt and hydroclimate-based runoff (P-E).
The amount of meltwater contribution reached up to 2 Mkm 3 per thousand years during the deglaciation (Fig. 10b). However, the timing and magnitude of meltwater contributions are highly sensitive to the isostatic adjustment (Figs. 7 and 8). When the isostatic component is included, the contribution of ice-melt was close to 58,302 m 3 s -1 at around 20 kyr BP, whereas, when isostatic component is not considered the ice-melt contribution reached a maximum of close to 36,263 m 3 s -1 at around 16 kyr BP (Fig. 10d).
Northward expansion of the Caspian Sea drainage basin, due to elevation change caused by ice-sheet growth, also resulted in considerable positive changes to P-E. The increase in P-E between 25 and 15 kyr BP is smaller when the isostatic component is considered (CSB-GIA) than when not considered (CSB-noGIA; Fig. 10c). The combined contribution from both ice-sheet meltwater and meteoric (P-E) water in the expanded glacial basin was substantial and likely the main contributor to the Caspian Sea water budget in the early deglaciation (Fig. 10d). The accumulated icemelt contribution was~17% higher than P-E from HadCM3 when the isostatic component was included (CSB-GIA), but~33% lower when its impacts were not included (CSB-noGIA). The treatment of the isostatic adjustment is therefore an important component to consider.
The P-E water budget we considered here uses the total amount of water falling as precipitation, and this may fall as rainfall or snowfall ( Fig. 10c and d). However, in reality, the climatic condition (very cold/dry or wet/warm) plays a significant role in how the precipitation is transformed to runoff. Snowfall may accumulate and retained as ice over prolonged periods and therefore not contribute to the water budget until the ice melts. Therefore, the climatic conditions over the ice-sheet play an important role in whether snowfall contributes to runoff. However, the model simulation output did not permit an analysis of this (only total precipitation was available), so we used the HadCM3 runoff field to integrate the various aspects of snowfall and melt (Fig. 10e). The modelled runoff reached its peak around~18e17 kyr BP and is 2000 m 3 s -1 less than P-E for CSB-noGIA (Fig. 10e). The amount of runoff from 25 to 19 kyr BP is roughly 50% less than P-E, likely as a result of a portion of the precipitation falling as snowfall that would not have contributed to runoff formation.
Runoff (HadCM3) and meltwater (ICE-G6_C) from the expanded basin contributed substantially to the deglacial Caspian Sea water budget (Fig. 10f). The ice-melt contribution dominates at around 16 kyr BP and 20 kyr BP for CSB-noGIA and CSB-GIA respectively, i.e. the timing is dependent on the isostatic adjustment. The accumulated ice-melt contribution (between 25 and 12 kyr BP) was higher than runoff from expanded Caspian Sea drainage area (by~42%) when isostatic adjustment was considered and lower (by~15%) when isostatic impacts were not considered.

Caspian Sea level variation using HadCM3 runoff, ICE-6G_C ice melt, and basin shape changes
The main drivers of Caspian Sea level variation during the last glacial period are thought to have been (1) climate change (P-E or runoff), (2) impact of ice-sheet growth/melt leading to restructuring of palaeo-drainage basin and reorganizing of the river systems, and (3) contributions of ice-sheet meltwater itself. In this section we integrate these different components into a simple model to calculate Caspian Sea Level in order to explore their relative importance (as described in section 3.2.4). We performed two sets of simulations, with and without bias-correction, using HadCM3 runoff, ICE-6G_C ice melt, and basin shape changes, identified as CSB-noGIA and CSB-GIA in Fig. 8. Results are presented in Fig. 11 and described below using HadCM3 runoff (output from land surface process model: MOSES2.1) to drive the model, but a similar figure is included in the supplementary information for comparison where HadCM3 P-E (output from the atmospheric component) is used instead of runoff (Fig. S4).
A first set of simulations were performed to evaluate the combined impact of runoff and ice-sheet meltwater in the expanded Caspian basin. As can be seen in Fig. 11a, the contribution of meltwater combined with runoff from the expanded area substantially increases Caspian Sea level starting from 25 kyr BP and lasting through the deglaciation period until 12 kyr BP. We set a sill height in the model at 40 m above m.s.l. (Manych-Kerch spillway). The increase in water budget resulted in overspill of the Caspian Sea to Black Sea basin between 20 and 15 kyr BP when the isostatic component is not considered (CSB-noGIA), and from 19 to 15 kyr BP when the isostatic component is considered (CSB-GIA). This can be compared to the simulation in which the basin is kept at the present-day extent (blue line in Fig. 11), where there is no large transgression leading to overspill at any point during the last 25 kyr. In this particular simulation, we note that although P-E (or runoff) over the whole basin is lower during the glacial period than during the Holocene (see Fig. 7a), the modelled Caspian Sea level is higher in the glacial period. This is because the model separates the water balance over the Caspian Sea and land, and there is a positive anomaly in P-E during the glacial period over the Caspian Sea grid points, which outweighs the lower runoff over land. We did a sensitivity analysis of the initial condition of the Caspian Sea volume (at 25 kyr BP) and found that it only impacted the Caspian Sea level evolution for~2 kyr, and did not have a long term influence (Fig. S5).
A second set of simulations were based on HadCM3 runoff integrated over the CSB-noGIA and CSB-GIA without ice-sheet meltwater contribution. As can be seen on Fig. 11b, runoff contribution from the expanded area alone greatly increased Caspian Sea level between 25 and 13 kyr BP when the isostatic component is not considered (CSB-noGIA), and between 25 and 16 kyr BP when the isostatic component is considered (CSB-GIA). The over-spill to the Black Sea occurred between 20 and 17 Kyr BP for CSB-noGIA, and from 19 to 17 kyr BP for CSB-GIA. The effect of the ice-melt is also to raise the height of the lowstand around 12 kyr BP by~23 m compared to when meltwater is not included (Fig. 11a and b).
In both sets of simulations, the impacts of meltwater input and runoff from basin expansion essentially decrease and become negligible by the early Holocene. In the Holocene all the scenarios are essentially the same, with a slight monotonic increase in Caspian Sea level throughout the Holocene, but much smaller variation than the deglaciation. The model reaches~28 m below m.s.l. of Caspian Sea level in the present-day even without bias correction. The simulated Caspian Sea level is~3 m less than the observed Caspian Sea level of~25 m below m.s.l. during pre-industrial period. This indicates that the modern-day model bias is small compared to the palaeo-changes. The bias correction (described in 13 section 4.2), which is applied over the whole period considered for the study has little impact on Holocene Caspian Sea level but does influence the magnitude of the extreme high and lowstands in the deglaciation, as shown in Fig. 11c and d.
In terms of the comparison between the model and palaeodata, both model (Fig. 11) and data (Fig. 6) suggest that there was a major transgressive episode leading to overspill into the Black Sea in the early deglaciation (19-14 kyr BP in the palaeodata). In the model, a key factor driving is the change in topography (due to ice-sheets) increasing the size of the basin. The increase in runoff/P-E into the Caspian alone is sufficient to lead to overspill in this time period (Fig. 11b). When ice-sheet meltwater is added this extends the timing of the connection to the Black Sea. The impact of isostatic adjustment opposes this and reduces the time interval of Caspian Sea overspill in the model (Fig. 11; compare dashed red line and black line). The impacts of Heinrich events (H1 in particular) on Caspian Sea level are more apparent when meltwater is not considered. Otherwise the impacts of meltwater on Caspian Sea overspill overwhelms the drying effect of the Heinrich events. The ice-sheet meltwater in the basin is immediately routed to the Caspian Sea and it is the derivative of the ice-sheet volume computed over thousand years with the assumption of a linear relationship (Fig. 10a). The actual variation in rate of ice-sheet melt during the thousand years was not considered as the ice-sheet reconstruction temporal resolution is very low. Therefore, what is shown in this result only reflects the volume of water that could have melted within the expanded area of Caspian Sea basin, leading to increased Caspian Sea level during cold dry periods (e.g. H1). Moreover, although the snapshot simulations (equilibrium climate simulations at every thousand years) allow us to effectively cope with constrained computational resources and boundary conditions, consideration of transient climate simulations would likely improve the last deglaciation period simulations characterized by abrupt climate change.
Both model and data indicate much more stable Caspian Sea levels during the Holocene, with the exception of the Mangyshlak regression in the early Holocene. The timing of the lowstand is earlier in the model, corresponding to the timing of the Younger Dryas, whereas in the data the average timing of the Mangyshlak is 9.5 kyr BP, although within uncertainties there is overlap. The magnitude of the lowstand recorded in the palaeodata during the LGM (Enotaevkan regression) is much lower than model output. In In the Holocene there is less variation in the model than in the palaeo-Caspian Sea level data. This is expected given the way the climate model was set up to simulate regular time slice intervals (which were run to equilibrium) that were processed as climatologies rather than a transient simulation with sub-millennial forcing (e.g. solar activity or volcanic variation) or internal variability on centennial/millennial time scales. This means we have not modelled some events that may well have produced Caspian Sea level variation at these time scales during the Holocene, such as the 8.2 kyr event (Alley et al., 1997) or the Roman climatic optimum, for example.

Discussion
In an attempt to identify the primary drivers of Caspian Sea level variation during the late Quaternary, we have investigated the nexus between hydro-climatological processes, ice-sheet evolution, and catchment dynamics and their potential impacts on Caspian Sea level variation. We integrated palaeoclimate model output from HadCM3 with reconstructed ice-sheet meltwater and topographic changes from ICE-6G_C to simulate Caspian Sea level over the last 25 kyr. We also collated Caspian Sea level data from published palaeo-environmental records for comparison.
Among the anticipated drivers, the Fennoscandian ice-sheet played the most substantial role on Caspian Sea level in the model by altering the topography of the earth surface and restructuring river drainage systems, increasing the area of the Caspian Sea drainage basin by 60e70% of its current size. The runoff in the expanded basin increased, despite the fact that the modelled glacial conditions were drier, due to the reorganization of river flows. When combined with meltwater during the last deglaciation period (19e12 kyr BP), these were the main drivers of the highest transgression of the Caspian Sea. The model only produced overspill into the Black Sea (as seen in the palaeodata) in the scenarios with an expanded basin area, and ice-sheet meltwater was not a necessary factor to produce the overspill in the model.
The climatic impacts of Heinrich events 1 and 2 and the Younger Dryas were included in the HadCM3 simulations. Each millennialscale event produces a decrease in the water budget (i.e. drier conditions over the basin) due to collapse of the AMOC. However, the impact of Heinrich event 1 on Caspian Sea level is overwhelmed by the influx of ice-sheet meltwater and runoff from the expanded Caspian Sea basin due to drainage reorganization at the time. This is not the case for the Younger Dryas, during which drier conditions lead to the lowest Caspian Sea level in the modelled time period. We note that, in the model, ice-sheet meltwater in the basin is immediately routed to the Caspian Sea. However, it is possible that the timing of meltwater routing may have been modified by proglacial lakes (Fig. 2[4]), which may have temporarily stored meltwater and then postponed and/or redirected its release (due to reaching overspill level, ice-damming collapse, or impact of isostatic adjustment). There is an indication in palaeorecords of pro-glacial lake water storage at the LGM and subsequent brief overflow into the Volga basin during the early deglaciation (Lysa et al., 2011;Mangerud et al., 2004). In addition, catastrophic flooding in Siberia (due to subglacial volcanoes) might have contributed to the overspill of proglacial lakes to the Aral Sea .
Moreover, as shown in Krinner et al. (2004), a number of large ice-dammed lakes, with a combined area twice that of the Caspian Sea, were formed in northern Eurasia during the last glacial period. Their result suggests that Eurasian proglacial lakes played a vital role in the reduction of ice sheet melting through strong regional summer cooling resulting in accelerated ice sheet growth and delay in ice sheet decay in Eurasia. Therefore, a careful consideration of the impacts of proglacial lakes on the Caspian Sea climate is necessary, and a spatially distributed hydrological routing model would be needed to explore the routing of the overspill that might have contributed to the Caspian Sea level.
In the literature, the timing of the Mangyshlak lowstand, and therefore the proposed cause, varies. In the model, the lowstand that resulted from Younger Dryas conditions overlaps in time with the uncertainties in the chronology of the Mangyshlak regression. The lowstand has been associated previously with the Younger Dryas (e.g., Bezrodnykh and Sorokin, 2016;Kislov, 2018) or the shift to warmer, dryer continental conditions in the early Holocene (e.g., Arslanov et al., 2016) or reduction in meltwater coupled with increased evaporation . HadCM3 (and indeed most other climate models) does not produce the magnitude of warming and increased continental conditions in the early Holocene that has previously been interpreted from the palaeo-record (termed the Holocene temperature conundrum by Liu et al. (2014). However, there has been recent debate about this interpretation of the palaeodata, following another study (Marsicek et al., 2018), which used fossil pollen temperature reconstructions over Europe and North America to demonstrate a gradual warming through the Holocene, with a similar temporal signature to climate models (including HadCM3), rather than an early Holocene climatic optimum. Such climatic conditions in the model are not able produce a significant lowstand in the early Holocene.
All of our hydroclimatic results are based on HadCM3, and our results could be sensitive to this choice. Kislov and Toropov (2007) show that a previous version of this model performed well compared to data for the Caspian for the present day and LGM. However, Bartlein et al. (2017) report that the ensemble mean of PMIP3 models show a drying of continental Eurasia during the mid-Holocene in contrary to data, which suggest a moistening. Examination of our simulations reveals that HadCM3 shows barely any change in precipitation over Europe (aside from some increase over Italy) during the mid-Holocene (Fig. S2), which is closer to the data. However, the data compilation has relatively few points over the Caspian Sea catchment and the available data shows some small drying, which is in better agreement with the model.
In the present day, there is no flow from the Amu-Syr Darya/ Himalayan drainage basins (Fig. 2) into the Caspian Sea, but there is geomorphological evidence that there have been previous shortlived connections that affected the water budget of Caspian Sea basin . Palaeo-runoff from the Amu-Syr Darya basin could have potentially contributed to past Caspian Sea level variations. However, as there is insufficient detail on the timing and extent of this connection, we were unable to fully assess any potential contribution. A brief analysis suggests that if the Amu-Syr Darya and the whole Himalayan basins flowed into the Caspian Sea, they could potentially increase the water budget by 2 and 5 times respectively (Fig. S1). However, contributions for the presentday leads to far too high P-E and doesn't reflect the current situation, and that available palaeodata suggests only brief periods of connection in the time period we've covered. Though there was potential for a large impact during periods of suspected connection, we think it is unlikely that the contributions were long-lived based on current evidence.
The Caspian Sea surface area changes rapidly with volume when Caspian Sea level is between À50 m and þ50 m above mean sea level. A larger surface area enhances evaporation from the Caspian Sea, adding more moisture in the atmosphere and affecting the water budget of the region as well as large-scale hydroclimate (Koriche et al., 2021b). We have made a first-order attempt to account for changes in evaporation due to surface area change of the lake in the simple model to simulate Caspian Sea level (see Equation (1)). In addition, the HadCM3 simulations have an increased lake size in the glacial compared to the present-day (see Fig. S2a), so there is an element of climate feedback included from the enhanced lake area, although this is prescribed and at the lowresolution of the climate model. Improved representation and coupling of all the processes that lead to changes in Caspian Sea area (e.g. ice-sheet meltwater) could further improve the Caspian Sea level model as well as the spatial patterns of change in the climate model.
As explained in section 2, there are competing ideas about sources and timing of runoff that caused the transgressive and regressive stages of Caspian Sea during the late Quaternary. One of the debates concerns whether the source of runoff that caused the Khvalynian transgression during the last deglaciation period was meltwater from the Fennoscandian ice-sheet via Volga river basin (e.g., Soulet et al., 2013;Tudryn et al., 2016) or glacial melting supplemented by proglacial lake outbursts from the Siberian region joining via the Aral Sea (e.g., Mangerud et al., 2004;Yanko-Hombach and Kislov, 2018). Based on our findings, the Fennoscandian ice-sheet played the vital role by altering the topography of the earth surface leading to river drainage system reorganization, and through meltwater contribution during the deglaciation period, although we cannot rule out some impacts related to proglacial lakes outburst and glacial melting from the Himalayan region. We find that the hydroclimate changes (P-E/runoff) at the LGM and early deglacial within the extent of the modern basin were not a dominant component in Caspian Sea level variation in our model when compared with ice-sheet related factors. This varies from the conclusions of some previous studies (Krijgsman et al., 2019;Yanchilina et al., 2019;Yanina et al., 2018). Holocene variations in Caspian Sea level were, however, dominated by hydroclimate change.

Conclusion
The Caspian Sea level variation is complex, resulting from different processes and settings. There are ongoing debates about the timing, mechanisms, and the source of runoff. Consequently, diverse thoughts exist on the importance of particular drivers of Caspian Sea level change. This study examined this further by evaluating the relative impacts of hydroclimatic change, ice-sheet accumulation and melt, and isostatic adjustment on the Caspian Sea level during the late Quaternary. Modelled Caspian Sea level variation, driven by the aforementioned forcings, was compared with newly collated palaeo lake level data covering the last glacial cycle. Overall, the results showed that the topographic change and the reorganization of the river drainage systems due to the Fennoscandian ice-sheet growth and retreat played the dominant role in the variation of the Caspian Sea level (especially the transgressive stages), although impacts related to pro-glacial lakes outburst and glacial melting from the Himalayan region may not be ruled out.
In conclusion the dominant forcing that controlled the Caspian Sea level during the deglaciation period (between 19 and 12 kyr) were the hydroclimate and ice-sheet accumulation/melting that resulted in the change in topography and glacial isostatic adjustment. However, during the Holocene the variations in Caspian Sea level were dominated by hydroclimate change. These results stress that source of the palaeo Caspian Sea level variations is complicated. Therefore, care should be taken when identifying the potential connections between the past and future sensitivities of the Caspian to climate change, as the glacial Caspian Sea level was partly influenced by activities that happened beyond its current catchment boundary and so are not relevant for the present-day and future. The interglacial periods are potentially more relevant than glacial periods for understanding the sensitivity of present and future Caspian Sea level variations. For example, the last interglacial period (130e115 kyr) climate (close to 4.0 C warmer than present with reduced ice-sheet extent globally) could be the best candidate to help understand the impacts of future climate change on the Caspian Sea. However, there are limited palaeoenvironmental records for validation and model initial (boundary) conditions for this time period. Putting aside those limitations, these results objectively identified the relative importance of different driving factors and their associated impacts on the extreme Caspian Sea level variations during the late Quaternary.

Data availability
The Climate model datasets used in this article can be found at http://www.bridge.bris.ac.uk/resources/simulations, hosted by Bristol Research Initiative for the Dynamic Global Environment of University of Bristol (Joy Singarayer and Paul Valdes, 2019).

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.