Inter-annual variability in snow cover depletion patterns and atmospheric circulation indices in the Upper Irtysh basin, Central Asia

,


Introduction
In the regions of Central Asia, owing to the arid and continental climate, snow cover melt provides an important contribution to runoff, determining the availability of freshwater for the local populations.Snow can, however, also act as a source of natural hazards, an issue which can be further exacerbated by climate change and human intervention (Dietz et al., 2014;Han et al., 2019).In Kazakhstan, snowmelt represents the primary source of water for agricultural and industrial use, dwarfing the contributions from liquid precipitation or groundwater; Kazakh rivers are characterized by large variability in annual runoff, with the largest annual discharges exceeding average annual runoff by 5-7 times and more than 75% of annual runoff occurring in a short period in early spring (Kozhakhmetov and Nikiforova, 2016).Floods caused by late winter and spring snow melt are one of the most severe natural hazards in the country, displacing 10-30,000 people each year and causing on average more than $30 million in damages (UNOCHA, 2016;Guha-Sapir et al., 2018) through loss of infrastructure, crops and livestock.Floods affect both lowland and mountain rivers, and there is evidence that the number of extreme events affecting mountain rivers has increased by 80% in recent years (Kozhakhmetov and Nikiforova, 2016).This issue is particularly relevant in the East Kazakhstan region, which was one of the regions most affected by flooding between 1967 and 1990 and the most affected between 1991 and 2015, with over 42 floods reported (Kozhakhmetov and Nikiforova, 2016).
In spite of the role of snowmelt for the hydrology of Kazakh rivers, comparatively little is known about large scale snow cover variations from year to year.While Mashtayeva et al. (2016) analysed the spatiotemporal distribution of snow depth and snow cover duration in all Kazakhstan, a detailed analysis of the patterns of snow cover depletion is necessary for the mitigation of floods, water management strategies and conservation of river ecosystems.The improvement of short-and long-term runoff models further requires an assessment of the relationship between snow cover variability, snowmelt and meteorological variables (Kuchment and Gelfan, 2007;Kang and Lee, 2014;Kult et al., 2014) and the prediction of timing and magnitude of peak runoff.Temperature exerts the greatest control on snow cover extent and duration (Bednorz, 2004;Hantel et al., 2000;Mote, 2006;Tang et al., 2017), although the spatial variation of the temperature-snow cover relationship with elevation is less well understood (Gurung et al., 2017).In addition, temperature anomalies and snow cover across Eurasia are related to large-scale atmospheric circulation patterns via feedback mechanisms (Cohen et al., 2012).Eurasian snow cover variability in different seasons has been linked with the phases of the North Atlantic Oscillation (NAO) and Arctic Oscillation (AO), Eurasian and Siberian pressure patterns and sea ice extent in the Barents-Kara Sea (Cohen & Barlow, 2005;Wegmann et al., 2015;Ye and Wu, 2017).Further still, spring snow retreat is also known to influence the East-Asian climate of the following summer by controlling the strength of the Indian summer monsoon (Zhang et al., 2017).However, the links between teleconnection patterns, temperature anomalies and snow cover retreat at the catchment scale and possible lagged effects of large-scale atmospheric circulation, which would provide further indications for long-range forecasts, remain to be fully explored.
The size of the East Kazakhstan region and the sparse network of meteorological stations (Mashtayeva et al., 2016) means investigations on snow cover cannot be realistically undertaken using field data, and leaves remote sensing as the most viable option.Information on snow cover extent from optical sensors such as MODIS (MODerate resolution Imaging Spectroradiometer) and AVHRR (Advanced Very-High-Resolution Radiometer) at moderate spatial resolution (500-1100 m per pixel) and at daily to weekly time scales has been employed for studying snow cover variability across Eurasia (Dietz et al., 2012;2013;2014).AVHRR data represent the longest temporal record from 1978 to the present, but requires additional processing for snow/cloud cover discrimination and the quality of the observations decreases significantly with the oldest generation of sensors (Dietz et al., 2014).Recently, the availability of Sentinel-2 data has opened up the possibility for investigating snow cover at a much higher spatial resolution (10-20 m) at weekly time scales (Hollstein et al., 2016), while Sentinel-3 continues the legacy of MODIS and AVHRR by providing data at moderate resolution (300 m) in large swaths.However, no long-term (i.e.decadal) records exist for the Sentinel satellites as the first was only launched in 2015.Alternative snow depth and snow water equivalent (SWE) datasets are available from passive microwave sensors (Pulliainen, 2006) or GRACE gravimetric data (Wang and Russell, 2016), though the spatial resolution of these products is much coarser, e.g. 25 km for the 'Globsnow' product (Luojus et al., 2013) and 1º for GRACE (Landerer & Swenson, 2012).In addition, passive microwave data are known to be unsuitable for snow cover detection in mountainous regions (Luojus et al., 2013), which form a large part of the study area.In this study, MODIS was chosen as the main data source, as it provides an easily accessible archive of snow cover data from 2000 to the present, produced through a robust methodology and requiring minimal additional processing (Hall & Riggs, 2007).
This study focuses on the analysis of snow cover variability in one of the main water catchments of Kazakhstan, using the MODIS MOD10A2 dataset.The aims of this study are: i) To investigate patterns of spring snow cover change in five sub-basins of the Upper Irtysh river catchment, including the magnitude and timing of early/late peak snow cover depletion rates and timing of snow cover disappearance at different elevations, and ii) to investigate the correlations between variability in snow cover, air temperature and atmospheric pressure patterns.The Irtysh basin is used as a case study to evaluate the general applicability of this approach to understanding the linkages between large scale snow cover change and runoff.
The Irtysh River, flowing in East Kazakhstan, is the main water resource of the country, with an average annual flow of 33.8 km 3 , and the most affected by floods (Kazhydromet, 2006).Its upper catchment (196,000 km 2 , Bayandinova et al., 2017;Fig.1)encompasses five sub-basins, namely Uba (9,917 km 2 ), Bukhtarma (13,663 km 2 ), Narym (1,852 km 2 ), Kurchum (4,975 km 2 ) and Kara Ertis (143,847 km 2 ) from northwest to southeast.Kara Ertis is the largest basin and the Kazakh name of the Irtysh river before the confluence into Lake Zaisan.The river originates in the southwestern Altay Mountains at 2500 m a.s.l and flows northwest through the uplands of northern China, where it is joined by streams originating in western Mongolia.The major tributaries of the upper river course originate on the western side of the Altay Mountains, and their basins are part of the East Kazakhstan administrative region.Narym, Kurchum (length: 218 km) and Bukhtarma (336 km) flow into the Bukhtarma reservoir; downstream of the city of Öskemen, the Irtysh is joined by the Uba River (Bayandinova et al., 2017).It then passes by Semey and Pavlodar (all with a population of 300,000 or above) in northern Kazakhstan and enters Russia, joining the Ob and eventually flowing into the Kara Sea.Water intake occurs in China and through a cascade of large water reservoirs namely the Ust-Kamenogorsk, Bukhtarma and Shulbinsk in Kazakhstan, (Huang, 2014), used for irrigation and hydroelectric power production; in fact, the transboundary nature of the river, shared between Kazakhstan, China, Mongolia and Russia, has led to conflicts over water usage (Hrkal et al., 2006;Alimkulov et al., 2017).
The climate of the study area is significantly affected by the presence of the Altay Mountains, capturing atmospheric moisture transported by the Westerlies and resulting in precipitation totals up to 1500-2000 mm yr -1 above 600 m a.s.l at the Irtysh headwaters (Zhang et al., 2018) and semi-arid rain-shadow conditions on the leeward Mongolian slopes (Huang, 2014;Klinge et al., 2003;Malygina et al., 2017).Temperature ranges are extreme, between +41°C in summer and -47 °C in winter (Malygina et al., 2017).In contrast with high precipitation in Upper Irtysh, the lower river basin is very dry, with annual rainfall not exceeding 250 mm at Pavlodar in North Kazakhstan, 500 km downstream of Öskemen (Huang, 2014).
3 Data Sources and Methods

MODIS Snow Cover
8-day composites of the MODIS MOD10A2 version 6 snow cover product (Hall & Riggs, 2016) were downloaded from the National Snow and Ice Data Centre.This product is derived from the temporal aggregation of MOD10A1 daily snow cover in 8-day composites and it provides the maximum snow cover extent (i.e. a pixel is classified as snow if snow is present on any day of the period).Snow classification is based on the 'snowmap' algorithm (Crane & Anderson, 1984), using the normalized difference snow index.MOD10A2 also includes information about lakes, lake ice and cloud cover.The first two are obtained using a lake mask, while cloud cover is reported if a pixel is cloudy in each daily image.The accuracy of snow classification in the MOD10A1 product has been reported as 93% (Hall & Riggs, 2007), and no further validation was performed in this study.We acquired tiles h22v03, h23v03, h23v04, h24v04 for each year between 2000 and 2017 and days of year (DOY) 57-201, corresponding to 26th February -20th July of a non-leap year, to enable analysis of the full snowmelt period (Table 1).Raw data tiles were reprojected from a MODIS sinusoidal grid to UTM45N WGS84, with uniform spacing of 500 m per pixel, mosaicked into a single raster layer and extracted for the area 40-54°N, 80-94°E.We processed yearly data using a cloud cover interpolation algorithm similar to previously reported in the literature (Gascoin et al., 2015;Tong et al., 2009) and calculated snow cover extents for each DOY and year (see Supplementary material).
Few data gaps were present in the original MOD10A2 and MOD10A1 datasets during 2001 (DOY 169 and 177) and 2002 (DOY 81), see also Dietz et al. (2013).Snow cover extent on these dates was calculated using linear interpolation based on the adjacent dates.

Elevation
We derived elevation information from the ASTER GDEM (Tachikawa et al., 2011), between 44-55°N and 80-94°E (a total of 165 individual tiles with 30 m spatial resolution).The DEM tiles were checked for inconsistencies and artefacts removed using a nearest neighbour approach.The tiles were then re-projected to UTM 45N WGS84, mosaicked to form a single large DEM and resampled to 500 m resolution to be consistent with the resolution of the processed MODIS snow cover product.

Basins
The basin shapefiles were rasterized and resampled to the same size as the 500 m MODIS grid cells.The rasters were then re-vectorized to provide areal information for a polygon corresponding to the grid size of MODIS data.Based on information from the ASTER GDEM, we extracted the hypsometry of each basin using 100 m elevation bins (see Fig. 2).

Calculation of snow cover variables
To investigate fluctuations and potential trends in snow cover in Upper Irtysh, we calculated the day of year (DOY) of snow cover disappearance (hereafter DSCD), which indicates when a pixel becomes snow free, and is similar to the concept of snow cover melt date used in other studies (Dietz et al., 2014;Liu & Zhang, 2017;Wang & Xie, 2009).We then averaged DSCD for each basin and, as snow cover is highly dependent on elevation (Dietz et al., 2012;Parajka et al., 2010), for each separate 500 m elevation range.
To provide further data for hydrological forecasting, we also calculated the peak snow cover depletion rate (hereafter PSCDR, expressed in km 2 d -1 ) by identifying the maximum snow cover difference between two images for each basin and elevation range and dividing it by 8 to calculate the depletion rate (equation 1).We similarly expressed the DOY of peak snow cover depletion (DPSCD) as the DOY when peak snow cover depletion occurs, according to equation ( 2).
Where SC is snow cover, h is elevation range, b is basin and the notation DOY + 8 and division by 8 in equation ( 1) is due to the 8-day composite period of MODIS.We also calculated values for each individual basin, by summing over elevation ranges.

Weather station data
We obtained air temperature observations from the National Oceanic and Atmospheric Administration (NO-AA) Global Historical Climatology Network (Menne et al., 2012).Three stations are located in the study area, all in Kara Ertis basin, namely Altay (737 m a.s.l.), Fuyun (827 m a.s.l.), and Baitag (1186 m a.s.l.).The stations are located south of the Altay Mountains following a northwest-southeast gradient (Fig. 1).Although daily minimum, maximum and average temperatures are available for these stations, we chose to use the latter as the maximum/minimum records are incomplete, often having less than 85% of available data each month.For the same reason, we discarded precipitation and snow depth observations.We chose 2000-2017 as our reference climate period to be consistent with MODIS data; daily mean temperatures were further averaged to obtain monthly, winter (DJF) and spring (MAM) records of mean air temperatures and calculate temperature anomalies relative to the reference period.

Reanalysis data
To further investigate the influence of atmospheric circulation on patterns of snow cover depletion in the whole study area, we downloaded ERA-Interim records of daily air temperatures and mean sea level pressure (MSLP) for central Asia from the European Centre for Meteorology and Forecasting, between 1990and 2017(Dee et al., 2011).ERA-Interim is produced by combining observational data with previous information from a forecast model in 12-hourly analysis cycles, and we chose it over other reanalysis products in view of its temporal coverage and relatively fine grid sampling of 0.75°x0.75°.We employed the daily analysis product with observations at 6:00:00 UTC, corresponding to noon local time, and averaged them to monthly, spring and winter values to calculate temperature anomalies over the reference climate period.

Atmospheric circulation indices
We obtained northern hemisphere monthly atmospheric circulation indices from the NOAA Climate Prediction Centre, which represent the leading patterns extracted from rotated empirical orthogonal functions (Barnston and Livezey, 1987), applied to 500 hPa height anomalies in the analysis region 20 -90°N.In turn, these reflect large-scale changes in atmospheric waves and jet stream patterns influencing temperature, precipitation and storm tracks as well as the position and intensity of the jet stream (NOAA, 2012).We included the indices that have been shown to influence the climate of Eurasia, specifically the NAO, the East Atlantic pattern (EA), the East Atlantic / Western Russia pattern (EAWR), the Scandinavia pattern (SCAND) and the Polar-Eurasian pattern (POL).
The NAO is the most prominent atmospheric mode during winter months, consisting in a north-south dipole of anomalies, one centred over Greenland and one over the North Atlantic between 35 and 40°N.It is associated with the intensity and location of the North Atlantic jet stream and storm track, with the positive (negative) phase leading to above-(below-) average temperature anomalies over Europe, extending over Siberia in prolonged phases (NOAA, 2012).The NAO has been widely used to analyse snow cover variability over Europe (Bednorz, 2004;Bojariu and Gimeno, 2003).The EA is the second most prominent mode, and consists of a north-south dipole of anomaly centres spanning the North Atlantic from east to west (NOAA, 2012).The positive phase of the EA is associated with above-average temperature in Europe and positive precipitation trends over Scandinavia for all months.The EA-WR pattern has four main anomaly centres: Europe, Northern China, north Atlantic and north of the Caspian Sea.It is a primary teleconnection affecting Eurasia throughout the year and in its positive phase it is commonly associated with below-average temperature anomalies for large portions of western Russia.(NOAA, 2012).The SCAND describes a primary circulation centre over Scandinavia and weaker centres over western Europe and western Mongolia.Its positive phase is characterized by anticyclonic anomalies over Scandinavia, giving rise to below-average temperatures over central Eurasia (Bueh and Nakamura, 2007).The POL is most commonly associated with above-average temperatures of eastern Siberia, reflecting an enhanced (weakened) polar vortex during a positive (negative) phase (NOAA, 2012).
We also included the Arctic Oscillation (AO), or Northern Hemisphere annular mode, which was obtained from NOAA CPC, based on 1000 hPa height anomalies.The AO is a planetary scale pattern of climate variability related to the zonal flow between 35 and 55°N, and has been widely linked to snow cover variability in Central Asia (Bamzai, 2003;Clark et al., 1999;Ye and Wu, 2017).Finally, we calculated the Siberian high index (SH), which expresses the strength of the Siberian high, a semi-permanent and quasi-stationary anticyclone stationed over Siberia throughout the winter and early spring months, associated with the coldest air masses of the northern hemisphere (Panagiotopoulos et al., 2005).The SH was calculated by averaging daily values of MSLP obtained from reanalysis data into monthly composites for the area 40-65°N, 80-120°E , as outlined by Panagiotopoulos et al. (2005) and standardized based on the 1990-2017 mean and standard deviation of MSLP for each month.For all indices, we produced winter and spring mean values by averaging DJF and MAM data.

Correlation analysis
We investigated links between snow cover variability, atmospheric circulation and temperatures through correlation analysis, assessing the strength of relationships using Pearson's correlation coefficient (r), and evaluating significance (p) using a two-tailed student test.We compared winter/spring indices with DSCD and DPSCD at separate elevation ranges and basins, and per pixel in the case of DSCD.To assess the influence of temperatures, we also compared winter/spring indices with temperature anomalies and temperature anomalies with snow cover variables.Temperatures from stations and ERA-Interim reanalysis between 2000 and 2017 were used for this test.Before this operation, data were detrended to isolate the inter-annual fluctuations (Panagiotopoulos et al., 2005).This analysis was also performed for each basin as a whole and at separate elevation ranges, with temperatures from each weather station and ERA-Interim grid cells in the study area.

Multiannual snow cover variability
While there are no notably strong trends in mean DSCD in any of the elevation intervals, there appear to be large inter-annual fluctuations at the lower elevations, especially in the larger Kara Ertis basin (Fig. 3d), while very little variability is seen above 3500 m a s.l..The most consistent patterns across the basins are the late DSCD events in 2010, best seen in Kara Ertis basin under 2000 m a.s.l., and in 2013, mostly seen in the smaller basins above 2000 m a.s.l.. Other notable features include the early DSCD in 2003, 2008 and 2012 above 2500 m a.s.l. in the smaller basins.2008 also shows a large shift towards earlier DSCD in Kara Ertis basin between 500 and 1000 m a.s.l..An isolated spike towards later DSCD is present in 2005 in Uba basin (see Fig. 3c), indicating snow did not completely melt during the observation period.At the lower elevations (< 1000 m a.s.l.), snow disappears on average before DOY 105 (mid-April, see Table 1), in the smaller basins and much earlier (DOY 73, third week of March) in Kara Ertis basin.At the upper elevations (3500-4000 m a.s.l.), it either disappears at the end of June in Kara Ertis basin, or does not melt at all during the observation period (until the end of July) in Bukhtarma basin (see Fig. 3b).Between 2500 and 3000 m a.s.l., snow disappears on average in May in Kara Ertis basin and early June in the smaller basins, with the exception of Uba basin in 2005.
In comparison with DSCD, larger inter-annual variability can be seen at all elevation ranges for DPSCD (Fig. 4).However, the patterns of DPSCD for entire basins seem to closely match those of their 500-1000 m elevation range, with the exception of Kurchum basin (Fig. 4a) where the basin-wide values lie between those of two elevation ranges, due to the different hypsometry of this basin (see Fig. 2).As with DSCD, late DPSCD can be seen in 2010 in Bukhtarma, Kara Ertis and Narym (see Fig. 4b,d,e) basins at most elevation ranges, while 2008 shows earlier DPSCD mainly at the lower elevations of Kara Ertis (< 2000 m a.s.l.) and Bukhtarma (< 1000 m a.s.l.) basins.The different elevation ranges generally show little crosscorrelation, with several contrasting values of early/late peak snow depletion, e.g. between 2500-3000 m and 3500-4000 m in Kara Ertis basin and 1500-2000 m and 2500-3000 m in Uba basin (see Fig. 4c).In the smaller basins, DPSCD generally occurs in April and on average, peak snow cover depletion occurs in early April in Bukhtarma and Uba (DOY 97) basins, and during the first week of April (DOY 89) in Narym basin.In the whole Upper Irtysh basin, DPSCD is synchronous with snow depletion in Kara Ertis basin except for 2010, when it occurred a week later.
PSCDR itself is largely controlled by snow retreat between 500 and 1500 m a.s.l., with smaller contributions from the other elevation ranges, and is up to 8 times larger in Kara Ertis basin, due to its greater size.For this reason, we here show the basin-wide values, not individual elevation ranges, together with the global value for the Upper Irtysh catchment (see Fig. 5).PSCDR for Upper Irtysh closely resembles the Kara Ertis basin series, with the exception of 2010, when there was a large drop in the Kara Ertis basin PSCDR.PSCDR maxima occurred in 2011 and 2017 with almost 5800 and over 5900 km 2 day -1 , respectively, with the lowest amount in 2007 (3160 km 2 day -1 ), and an average of over 4500 km 2 day -1 .We evaluated trends for each basin and the whole Upper Irtysh using both ordinary least-squares linear regression and the Mann-Kendall test.A significant (p < 0.05) positive trend of +67 km 2 day -1 per year was identified for the Upper Irtysh basin.However, the trend is only significant at the 90% confidence level when the Mann-Kendall test is used.No other significant trend is evident for individual basins.

Atmospheric circulation indices and temperature anomalies
To identify the most suitable indices for snow cover analysis, we first looked at their correlation with temperatures, by considering aggregated winter (DJF) and spring (MAM) indices.The AO and SH emerge as the strongest indices: both the winter AO and SH are correlated with temperature anomalies from all three stations and across basins in the reanalysis data, although the AO is non-significantly (p > 0.05) correlated at some grid points (Fig. 6a).Specifically, the winter AO appears positively correlated with winter temperatures across Siberia exerting its influence north to the Barents Sea and south to the Altay Mountains (see Fig. 6a).Conversely, the Siberian high exerts a more local impact in winter, influencing temperatures north of the Tibetan plateau, across Northern China, East Kazakhstan and Mongolia, with its southernmost tip over southeast China (Fig. 6c).
Spring SH and winter AO are also strongly correlated with spring temperature anomalies at all stations and grid cells in the Upper Irtysh (see Table 2 and Fig. 6b,d).The correlation between the winter AO and ERA-Interim spring temperature anomalies is significant (p<0.05)across the Altay Mountains, between eastern Kazakhstan, northern Mongolia and southern Russia, including the study area (see Fig. 7b), with patterns of significant correlation with temperatures displaced to the east compared to the winter.The correlation between spring SH and spring temperature anomalies from ERA-Interim shows a similar pattern to the correlation between winter temperature and winter SH, although the extent of the influence is limited to a latitudinal band between 70°E and 120 °E (see Fig. 6b,d).
Other indices show moderate correlation at either individual stations or grid cells, e.g. the winter NAO is also correlated with winter and spring temperatures, although the strength of the relationship is generally lower than for the AO, except for the winter index -winter temperature anomalies at ERA grid cells.The winter SCAND and EA also correlate with winter temperatures.The first correlation however is only significant at weather stations, while the latter at two stations and limited grid cells (see Table 2).

Atmospheric circulation indices and snow cover variability
The relationship between the two strongest indices (i.e.winter AO and spring SH), spring temperature anomalies from reanalysis data and Altay weather station and DPSCD in the Upper Irtysh are plotted in Fig. 7.All variables show significant moderate to strong cross-correlation, particularly temperature anomalies from ERA-Interim and Altay stations (0.94).Additionally, high interannual variability exists for the AO and temperatures, while the SH is mostly negative from 2000 to 2017.The largest anomalies in all series occur in 2010, with negative values for the AO and temperatures and positive ones for the SH and DPSCD.In contrast, positive anomalies in temperatures in 2008 are synchronous with a positive AO and negative SH and DPSCD.In spite of the general agreement between the variables, some exceptions occur, such as positive temperature anomalies in 2004 followed by a late DPSCD and the negative winter AO in 2013 with above average temperatures the following spring.
Among the two indices, the winter AO shows the strongest (negative) correlation with DPSCD in all basins, while Spring SH is significantly positively correlated with DSPCD only in Kara Erys and Narym basins.Both indices also correlate with Upper Irtysh basin as a whole (see Fig. 8 and Table 3).Across different elevation ranges, winter AO shows moderate to strong correlation with DPSCD up to 1500 m a.s.l. in Kara Ertis, Bukhtarma and Narym basins, with generally higher values in the Kara Ertis basin (max |r| 0.73) compared to the others (see Fig. 8a,c).As regards spring SH, the highest (positive) correlations are again seen in Kara Ertis basin up to 2000 m a.s.l..In the other basins, correlation is either not significant or significant in fewer elevation ranges compared to the AO (see Fig. 8b,d).
The per-pixel correlation between the AO and SH and DSCD shows stronger correlations compared to DPSCD across multiple elevation ranges (see Fig. 9).Spatially, the correlation with the winter AO shows elevation dependence, with the strongest values below 2000 m a.s.l.south of the Altay Mountains in most basins.The main exception is an area of lower correlation in the south-eastern part of Kara Ertis, on the Mongolian side, corresponding to the dry steppes of the Gobi Desert (see Fig. 9a).Across elevation ranges, moderate to strong significant negative correlation with winter AO is seen up to 2000 m a.s.l.(3500 m a.s.l. in Kara Ertis basin), with the strongest values all between 1000 and 1500 m a.s.l.(see Fig. 9c).The spring SH shows generally lower correlation values, mainly limited to the Kara Ertis basin, while very few pixels correlate in the smaller basins, limited to the lower elevations.The highest correlations are seen in the Gobi Desert, the same area that showed lower correlation with the winter AO, and west of lake Ulungur in the southern Kara Ertis basin (see Fig. 9b).Moderate significant positive correlation with DSCD is seen up to 2000 m a.s.l. in this basin.Significant correlation is also seen up to 1000 m a.s.l. in Narym and 500 m a.s.l. in Bukhtarma and Uba basins (see Fig. 9d).
Both snow cover variables show significant correlations with spring temperature anomalies from weather stations in the relevant elevation band, i.e. 500-1000 m a.s.l. for Altay and Fuyun and 1000-1500 m a.s.l. for Baitag.For DPSCD, the strongest correlation is seen in Kara Ertis basin with a max.|r| of 0.79 at Altay; this basin has the highest correlations at all three stations, while no significant correlation with temperatures from any weather station is seen in Uba basin (see Table 4).For DSCD, all correlations with temperatures from weather stations are statistically significant: the strongest correlation is between DSCD in Kara Ertis basin and Baitag spring temperature anomalies, with high values also in Narym and Uba basins (see Table 4).

Data and methods uncertainty
The main uncertainties in the analysis of snow cover variability and the association between snow cover variables and atmospheric indices are related to 1) the use of 8-day snow cover area instead of daily products, which also implies the derived DSCD and DPSCD are in fact 8-day products, and 2) cloud-cover interpolation.Drawing links between these variables and flood hazards in East Kazakhstan also requires an assumption that snow covered area can be used to represent SWE and snow cover depletion can be a proxy for snow cover melt and river discharge.
The choice of the 8-day composite product was made to reduce the impact of cloud cover on the study area.While the main misclassification in the daily MODIS product occurs between clouds and snow, and can propagate to the 8-day product (Hall & Riggs, 2007), the accuracy of the composite product might actually be higher than that of the daily products (93% on average, Hall & Riggs, 2007) in areas of high cloud cover, also leading to a higher correlation of the composite product with streamflow (Zhou et al., 2005).In a study conducted in northern China, Wang et al. (2009) found an accuracy of 94% for snow mapping for the composite product at snow depths greater than 4 cm, decreasing to 39% for lower depths and patchy snow cover.These are however likely to represent a small amount of SWE and therefore have a small impact on the assessment of water availability.The use of a cloud-interpolation scheme can also help reduce the cloud-cover related uncertainty (Dietz et al., 2011;Dong & Menzel, 2016), although our algorithm might introduce biases by enhancing elevation dependency (see S1, point 3) and in case of snow melt and subsequent snowfall (see S1, point 4); a further improvement on cloud cover interpolation accuracy is the use of in situ temperature and precipitation data (Dong & Menzel, 2016), although these are not often readily available in remote areas.Small uncertainties in snow cover variability might also stem from data gap filling in the original MODIS product in 2001 and 2002 at lower and higher elevations, respectively.Higher elevations in our study area might be affected by late spring snow falls after snow melt, which will not be recorded as DSCD, calculated as the first week of snow disappearance.However, these events are expected to provide a low contribution to the SWE.
The association between snow covered area, SWE and stream discharge has been demonstrated in several studies using the daily or 8-day MODIS composite product and other SCA datasets.Delbart et al. (2015) found an uncertainty of 15% in predicting water discharge from SCA between April and September in four Andean catchments.Further still, Gurung et al. (2017) found high correlations (> 0.78) at elevations between 4000 and 6000 m a.s.l. in the Gandaki basin in the Hindu Kush Himalayas.Using the composite product, Tong et al. (2009) identified a correlation of 0.84 between cloud-filtered snow cover extent and streamflow in the Quesnel river basin, Canada.
In this study, the relationship between snow covered area and SWE cannot be assessed due to the lack of in situ observations in our study area (Mashtayeva et al., 2016).The regulated nature of the Irtysh River, i.e. the presence of reservoirs and channels in China and Kazakhstan also causes dampening of peak flow and decreased runoff variability downstream (Huang et al., 2012;2014), preventing attempts at correlating snow covered depletion rates and timing and discharge.

Variability and trends in snow cover depletion
In this study, we found no significant trend but strong fluctuations in DSCD and DPSCD in all sub-basins of the Upper Irtysh between 2000 and 2017.Both variables appear strongly negatively associated with spring temperature anomalies at elevations below 1500 m a.s.l.(see Table 4), which represent the largest contributing areas to runoff for most basins (see Table 2), with the highest correlation in the largest sub-basin of Kara Ertis.Above these altitudes, temperature variations are likely less relevant, and other topographic features such as aspect.shading, snow redistribution by avalanching or wind may dominate over the influences of regional climatic variability.Fluctuations in DSCD agree with the findings of Dietz et al. (2011) in their study of snow cover in central Asia using MODIS data, with matching years of relevant snow cover disappearance anomalies (e.g. 2008, 2010).High multiannual variability in Kazakhstan snow cover during the melting season was also found by Mashtayeva et al. (2016) although snow depth was the variable investigated.Considering longer periods of observation, Dai & Che (2014) found a significant decrease in late spring snow depth and snow cover duration (-1 day/year -1 ) in the Chinese Altay between 1987 and 2011, while in the Ob basin, Yang et al. (2003) identified an upward trend in spring snow cover between 1966 and 1999, caused by increased precipitation over Western Siberia.Long-term temperature records between 1966 and 2015 from southern Altay show warming trends in winter and spring, with the greatest rise in spring (1.15°C per decade at Fuyun, Zhang et al., 2018).Warming in average annual temperatures were also reported by Hu et al. (2014) using several reanalysis datasets between 1979 and 2011.If the warming trend continues as projected, a shift towards earlier snow cover disappearance and peak snow cover depletion would be expected to occur in the Upper Irtysh and across Eurasia.The existence of longer-term trends for the study area might be investigated using AVHRR data, albeit with lower spatial accuracy (Dietz et al., 2014).
In this study, peak snow cover depletion rates were not significantly correlated with the timing (DPSCD) in any sub-basin.PSCDR are also uncorrelated with the onset of melting conditions at weather stations, including rates and magnitudes of temperature increase at crossing of the 0º C isotherm.A significant correlation (0.50) was found between PSCDR and snow cover area on DOY 57 in late winter (see Table 1), which can represent a proxy for maximum SCA for a given year.Late winter snow cover explains 25% of the variance in PSCDR, but it shows no significant trends between 2000 and 2017, thus it does not appear to be the cause of the observed increase in PSCDR (Fig. 5), which remains unresolved.Among the sub-basins, Narym and Kurchum show the highest multi-annual variability; large peaks are also seen in Bukhtarma (2001 and2010) and Uba (2000 and2011) basins, with potential implications for local flood risks and water availability for reservoir operation.

Atmospheric circulation indices and snow cover depletion
Our correlation analysis shows the AO and SH to be the circulation patterns most influential to patterns of snow cover depletion in the Upper Irtysh basin in winter and spring.The correlations between winter AO anomalies and DPSCD are significant below 1500-2000 m a.s.l.: notably, DPSCD shows higher correlation with the winter AO at higher elevation, while the opposite is true for the SH, and a similar inverse pattern is seen for DSCD.
Two winters stand out as highly anomalous in both the AO and SH record during the period of observation: 2007-2008 (positive) and 2009-2010 (negative), with matching winter temperature anomalies and patterns of subsequent spring snow cover depletion (Fig. 7).Winter 2009-2010 was in fact remarkably cold across Europe and central Asia.According to Cohen et al. (2010), the highly negative AO was initiated by the rapid advance of Siberian snow cover in October, followed by the strengthening of the SH, increased upward wave activity flux, stratospheric warming and downward propagation of height and wind anomalies to the surface.The link between autumn snow cover and phases of the AO has been reported in several studies, with feedback mechanisms including soil moisture and the albedo (Bojariu & Gimeno, 2003, Cohen & Entekhabi, 1999, Saito & Cohen, 2003); patterns of decline in Arctic sea ice have been linked to the AO (Liu et al., 2012).Compared to 2009-2010, the winter of 2007-2008 received little attention, although unusually mild temperatures were reported across Eurasia, linked with the positive phase of the NAO/AO (Ilkka et al., 2012), in addition to a La Niña event (Stevenson, 2016).
In our study, the winter AO AO is strongly negatively correlated with DPSCD and appears to influence spring snow cover retreat through the modulation of spring temperature anomalies.The lagged influence between the AO and spring snow cover has in fact been previously reported in the literature.Bamzai (2003) observed a correlation of 0.46 between JFM AO and AM snow cover melt date over Eurasia.Foster et al. (2013) found that the winter AO is negatively correlated with spring snowmelt and DSCD at latitudes between 50º and 60º N, while the opposite is true between 60º and 70º N. In Saito & Cohen (2003), evidence for reciprocal lagged effects was found, i.e. late winter/spring AO leading spring/summer snow and spring/ summer snow leading the following fall/early winter AO.In a subsequent study, Saito et al. (2004), found a consistent relationship between fall Eurasian snow cover and the winter AO between 1971 and 2000.During the 1980s, an association between winter AO and spring snow cover also emerged, as well as autocorrelation between winter and spring snow cover and winter/spring AO.This latter correlation was also reported by Foster et al. (2013) but found no confirmation in our study.While our analysis did not include fall snow cover, its possible direct influence on winter AO and spring snow cover retreat is appealing for long-range forecasts and should be discussed in future research, both in the Upper Irtysh and other Eurasian catchments.
Compared to the NAO, the winter AO exhibits a stronger connection with both winter/spring temperature anomalies and DPSCR/DSCD.While the similarity/difference between the AO and the NAO and the actual physical meaning of the AO are still a matter of debate (Báez et al., 2013), our study confirms the hypothesis of a weaker (Saito & Cohen, 2003) and more local effect of the NAO on snow cover variability, restricted to Europe and western Asia (Clark et al., 1999).The observed association between the spring SH, temperature anomalies and DPSCD, while confirming that the Siberian high is indeed a relevant atmospheric influence on the climate of Kazakhstan (Panagiotopoulos et al., 2005), also show that it holds a limited potential in hydrological forecasting, as no relevant lagged correlations were found in our study.In spite of evidence that a positive SH reinforces negative AO anomalies (Cohen & Entekhabi, 1999, Cohen et al., 2010), the interplay between the AO and the SH is not entirely understood (Huang et al., 2016), and whether one index can be predicted from the other still needs to be assessed.

MODIS and the AO in flood forecasting
High interannual variability in timing and rate of peak snow cover depletion appear to modulate flooding occurrence, intensity and timing in eastern Kazakhstan.Specific flood reports for the study area are lacking, but figures from the Emergency Events Database (Guha-Sapir et al., 2018) and Global Active Archive of Large Flood Events (Dartmouth Flood Observatory, 2018) point to snowmelt-induced floods occurring throughout Kazakhstan in 2008, 2010, 2011, 2015and 2017 (Fig.5).Comparison with peak PSCDR shows a general agreement, with the highest peak snowmelt rates in the 18 year series occurring in 2011 and 2017.In particular, severe flooding affecting the Uba basin in mid-April 2011 was reported on the official government news website (Dixinews, 2011), corresponding with the timing of the highest PSCDR in the Uba basin in the study period (Figs. 4 and 5).However, the third highest PSCDR in 2006 has no correspondence in flooding, although this might be due to a bias towards better reporting in more recent years (both databases are based on local news).Broad agreement is also seen with respect to timing, with floods occurring relatively late in 2011, 2015 and 2017 (end of March -beginning of spring) and early in 2008 (20th February, Dartmouth Flood Observatory, 2018), corresponding with late and early values of DSCD (Fig. 7).Flood occurrence in 2010 however is reported as early as March, in spite of very late peak snow cover depletion found in our study.2010 was a particularly cold winter with record snowfalls (Cohen et al, 2010), and temperature increase above the snow melting point at low elevations might have led to abundant snow melt and subsequent flooding.
In spite of this general agreement, quantitative assessment of the predictive ability of snow cover depletion maps based on MODIS and of the possible use of the AO in flood forecasting in the Upper Irtysh catchment is complicated by a lack of in situ observations in the study area.To improve the ground-based hydrological monitoring network in Kazakhstan, runoff should ideally be measured at stations, located in the smaller, nontransboundary basins upstream of Bukhtarma reservoir.In principle, the timing of flooding could could also be investigated via remote sensing, by detecting flooded areas through optical (Revilla-Romero et al., 2015) or SAR (Brown et al., 2016) satellite images, and algorithms have also been devised to estimate discharge from MODIS images based on the reflectance difference between water and land pixels in near infrared bands (Tarpanelli et al., 2017).Validation of these approaches in sparsely gauged catchments however remains a non trivial task.Since our approach is based on readily available data, in areas where ground based observations do exist, it could be more easily validated, providing a means to evaluate MODIS DPSCD and PSCDR products with respect to water management strategies and possibly implement long-term forecasts based on the AO.

Conclusion
In this study, we explored multi-annual (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017) snow cover variability in the upper course of the Irtysh, a large transboundary river basin of Central Asia, which is significantly impacted by spring snowmelt floods.Given the scarcity of direct runoff observations in the region, we analysed snow covered area as a means to understand the seasonal and inter-annual variability of snow hydrology, to provide valuable information in the context of flood forecasting.The analysis was carried out using 8-day MODIS snow cover composites, both within the five major tributary basins and the Upper Irtysh as a whole.Cloud filtered data were used to investigate: 1) the timing of snow cover disappearance, using the day of year when snow cover disappears from different 500 m elevation ranges within each basin (DSCD), 2) peak snow cover depletion, the day of year of maximum areal snow cover depletion in a basin (DPSCD) and 3) peak snow cover depletion rate in each basin in each year (PSCDR expressed in km 2 d -1 ).In addition, we employed ERA-Interim reanalysis data, local weather stations and atmospheric circulation indices to investigate potential meteorological drivers of the observed snow cover depletion patterns.
The analysis of snow cover variability in the Upper Irtysh basin points to large inter-annual and inter-basin differences, for both DSCD and DPSCD, with the largest, low-lying basin showing the earliest dates and providing the largest snowmelt contribution to runoff.DPSCD generally occurs between mid-March and April, depending on the basin, while DSCD occurs later during spring.No clear trends over the 2000-2017 study period were identified in any basin for either variable.In contrast, PSCDR, which peaked at over 5900 km 2 day -1 in 2017, displays a weak increasing trend for the Upper Irtysh basin over the same period.In view of their relevance for flood hazards, the existence of longer-term trends, as well as the causes of the increase in PSCDR, should be addressed in further research.
Multi-year fluctuations in peak snow cover depletion and snow cover disappearance are controlled by spring temperature anomalies and regional atmospheric patterns detected by weather stations and reanalysis data in the Upper Irtysh basin.The atmospheric patterns that best correlate with temperatures and snow cover variables are the winter Arctic Oscillation and spring Siberian High indices.Among the two, the winter Arctic Oscillation shows the strongest correlation with DPSCD, especially in the larger, low-lying basins (-0.71), while lower correlations are seen in the smaller mountain basins, where topographic variability likely plays an important role.
The ability of the winter Arctic Oscillation index to predict a relatively late or early timing of the spring peak snow cover depletion rate in parts of Eurasia makes it a potentially useful tool for long-term forecasts, and for use in early-warning flood prevention schemes.However, the lack of empirical runoff data in the Upper Irtysh basin complicates validation of the feasibility and effectiveness of this approach.Despite the lack of validation, our approach is based on readily available data and can therefore be applied to similar catchments where snowmelt is the dominant source of runoff, providing meaningful insights into the interactions between large scale snow cover depletion patterns and streamflow, for informing water management strategies.A denser network of ground based observations or different remote sensing approaches could be employed to further assess the influence of the Arctic Oscillation in modulating snow cover depletion at the catchment scale.

Acknowledgments
This research was supported by the Newton Al-Farabi Industry-Academia Partnership Programme (project number IAPP\1516\42) administered by the Royal Academy of Engineering (UK) and a Student Mobility

Table 1 :
Day Of Year (DOY) -date correspondence for 8-day composite MOD10A2 data.Dates are shown for a non-leap year.

Table 2 :
correlations between average winter (DJF) and spring (MAM) indices and DJF -MAM temperature anomalies from weather station and ERA-Interim data between 2000 and 2017.Bold italics indicate significant correlation (p < 0.05) at individual stations/ERA-Interim cells and indices with highest correlations.