Changing snow seasonality in the highlands of Kyrgyzstan

Snow cover extent has been observed to be changing for more than three decades using both in-situ data (Groisman et al 2006, Bulygina et al 2010, 2011) and remote sensing products (Schanda et al 1983, Hall et al 2002, Brown and Robinson 2011). Remote sensing techniques for the monitoring of snow cover extent have advanced substantially in the past three decades (Robinson et al 1993, Armstrong and Brodzik 2001, Painter et al 2009, Rittger et al 2013, Morriss et al 2016). Much of the change analysis has focused on broad scales—from hemispheric to subcontinental—at coarse spatial resolution (Robinson and Dewey 1990, Groisman et al 1994b, Brown 2000, Robinson and Frei 2000, Dye 2002, Déry and Brown 2007, Hori et al 2017).

The impacts of climate change are exacerbated in mountainous regions (Beniston 2003, Fischlin et al 2007, Immerzeel et al 2010, Rangwala and Miller 2012, Kohler et al 2014). However, the ability of global climate projections to capture the complex dynamics of mountain climates is limited (Christensen et al 2013), particularly over mountainous Central Asia (Hijioka et al 2014, Reyer et al 2017).

Studies in montane Central Asia are relatively few, especially at higher spatial resolutions, but they have all shown significant changes, whether in snow cover (Dietz et al 2013, 2014, Zhou et al 2013, Tang et al 2017), glacial extent (Aizen et al 1995, Narama et al 2010), or meltwater runoff (Aizen et al 1997, Chevallier et al 2014). Still, there is a notable paucity of studies on the changing environment across montane Central Asia (Hijioka et al 2014, Reyer et al 2017).

Here we evaluate trends in the seasonality of snow cover across Kyrgyzstan since 2002 using the latest version of the MODIS (MODerate Resolution Imaging Spectrometer) snow cover composites. Our need to quantify changing snow cover seasonality arises from our interest in highland pasture conditions in rural Kyrgyzstan. Kyrgyz livelihoods based on montane agropastoralism are particularly vulnerable to changing environmental conditions due to a reliance on the seasonal movement of livestock to higher elevation pastures, a practice also called vertical transhumance (Schillhorn van Veen 1995). We assess at multiple scales whether and where the snow cover timing and duration have changed significantly in the recent past: at the scale of oblasts (provinces); in four rayons (districts); and at elevational bands within these selected rayons. We compare results from MODIS snow cover products generated from different overpass times (late morning and early afternoon), and two product versions.

Zoom In Zoom Out Reset image size

The study area is the territory of Kyrgyzstan (Kyrgyz Republic), a highly mountainous country in the eastern part of Central Asia. Kyrgyzstan shares borders with, moving clockwise from due north: Kazakhstan, China, Tajikistan, and Uzbekistan. The total area of the country is just shy of 200 000 km², of which 191 801 km² is inland and 8150 km² is in open water. In 2017, the population was approximately 5.8 million living in seven oblasts or provinces (figure 1(a)), of which only 36% are in urban areas with the largest city being the capital Bishkek.

Kyrgyzstan is among the poorer nations, with an estimated per capita gross domestic production based on purchasing power parity of just US$3551 in 2016 (WorldBank 2017a). Moreover, the country is heavily dependent on remittances from workers aboard; the World Bank projects remittance inflows of US$2.5 billion for 2017 (WorldBank 2017b). Poverty, particularly rural poverty, limits adaptive capacity to respond to impacts on livelihoods arising with climate change (Lioubimtseva and Henebry 2009, Reyer et al 2017).

More than 56% of the territory occurs above 2500 m (Azykova 2002). Mountain ranges cover more than 90% of the land area. These ranges include parts of the Pamir and the Alatau, and a large portion of the Tien Shan that divides the country into two zones. The northern zone holds three oblasts—Talas, Chuy (including the capital Bishkek), and Issyk-Kul—and the southern zone has four oblasts—Jalal-Abad, Naryn, Osh, and Batken (CACILM/ADB 2010).

The climate of Kyrgyzstan is influenced by country's inland location between temperate and sub-tropical zone, high elevation, the distance from oceans, and proximity to the deserts. It results in intense solar radiation, low precipitation, and a continental climate (Akimaliev et al 2013). The mountain relief causes elevational climatic zonation of temperature and moisture. In the hot months of July and August, the mean air temperature over lowlands ranges between 17 °C–40 °C and only ~4 °C in the mountains. During winter months, the lowest temperatures are recorded in the mountain valleys and depressions (Kulikov and Schickhoff 2017), but frost can occur in every oblast. Annual precipitation varies from 144 mm in some parts of Issyk-Kul to 1090 mm in the lowlands of the Fergana valley, but precipitation is unevenly distributed across the country. Vegetation types are scattered along distinct elevational zones, influenced by vertical gradients of climatic variables. Less than 10% of land area is appropriate for crops, forests cover ~5%, and more than 50% of the land is used as pastoral rangeland (CACILM/ADB 2010).

3.1. Satellite data

To characterize snow seasonality, we used the most recent Version 6 of the MODIS Terra snow cover 8 day composites with a nominal spatial resolution of 500 m (MOD10A2/MYD10A2) for the period of 14 years starting in 2002. Both datasets (Terra and Aqua) report the maximum snow cover extent observed during an 8 day period by compositing 500 m observations from the MODIS daily snow cover products (MOD10A1/MYD10A1) generated by The National Snow and Ice Data Center (https://nsidc.org/). Snow cover information is derived from the Normalized Difference Snow Index (NDSI) (Hall et al 2002). Snow cover typically has very high visible (VIS) reflectance and very low shortwave infrared (SWIR) reflectance. The algorithm uses a threshold test for spectral band ratios of a difference in VIS (band 4: 0.555 μm) and SWIR (band 6: 1.650 μm) reflectance.

$$\begin{equation} {\rm{NDSI}} = \frac{{\left( {{R_{{\rm{vis}}}} - {R_{{\rm{SWIR}}}}} \right)}}{{\left( {{R_{{\rm{vis}}}} + {R_{{\rm{SWIR}}}}} \right)}}. \end{equation} \tag{ 1 }$$

NDSI > 0.0 indicates the presences of some snow within the pixel, while a pixel with NDSI < 0.0 indicates a snow-free land surface (Riggs and Hall 2015).

Although snow cover always has NDSI > 0.0, not all surface features with positive NDSI values are snow covered (i.e. salt pans or cloud-contaminated pixels at cloud edges). Thus, additional screening procedures are required to reduce commission error. In the MODIS products, a pixel will be mapped as snow if the NDSI > 0.4 and the MODIS band 2 (red) reflectance exceeds 0.11 to discriminate snow from water (Riggs and Hall 2015). These MODIS snow cover products (MOD10A1/MYD10A1) have undergone extensive peer-review evaluation and validation processes (Hall et al 2002, Riggs and Hall 2015), and no effort was made in this retrospective study to conduct additional product evaluations.

We postprocessed the products by mosaicking two MODIS tiles (h23v04 and h23v05), reprojecting into WGS-1984, and extracting all pixels flagged as 'snow'. There were 46 composited images per year per product. We chose to work with the 8 day composited products instead of the daily products for two reasons: (1) so that the statistical power for the trend analyses would be constant across scales, and (2) because far fewer studies have worked with the composited data.

For terrain information, we used the 15 arc-second (~450 m) mean elevation product of Global Multi-resolution Terrain Elevation Data 2010 (GMTED2010) developed by the US Geological Survey and the National Geospatial-Intelligence Agency.

3.2. Trend analysis

To analyze changes in snow cover seasonality, we defined the observation season for each year to start on day of the year (DOY) 169—approximately the summer solstice—and extend to DOY 168 of the following year (DOY169_year through DOY168_year+1). Thus, we analyzed a 14 year time series starting in the middle of 2002 and ending in middle of 2016. For each snow season, we tracked three snow cover variables: the first date of snow (FDoS), the last date of snow (LDoS), and the duration of snow season (DoSS). FDoS was the composite date when the snow pixel is marked as 1 (snow on) first time for each snow season. LDoS was calculated inversely, the composite date of the last appearance of snow during the snow season. The DoSS was the simple difference between the LDoS and the FDoS. For each of the snow variables, we calculated the mean, standard deviation, and coefficient of variation for the series of 14 snow seasons from 2002/03 through 2015/16.

Note that since we are interested in the potential impact of changes in snow seasonality on pastoralism, we have purposefully chosen to characterize the snow cover season by its temporal extremities: the first occurrence of snow appearing in a composite during the observation season (FDoS) and the last occurrence of snow appearing in a composite (LDoS). We understand that snowmelt can occur after FDoS and before LDoS, even multiple times. Were our purpose to evaluate snow cover duration to estimate the regional hydrological budget or the surface energy balance, then these outer bounds of snow occurrence could overestimate snow cover influence. However, our motivation in characterizing snow season timing is different. We are more interested in pasture dynamics than in high mountain snow processes.

Simple linear regression has been used by remote sensing scientists to estimate trends (de Beurs and Henebry 2008). However, it is better to use non-parametric tests since they provide higher statistical power in case of nonnormality and are robust against outliers (de Beurs and Henebry 2004). To evaluate the change in snow season metrics, we applied the non-parametric Mann–Kendall trend test and the Theil-Sen linear trend estimator. The non-parametric test is based on the rank correlation coefficient statistic τ (Kendall 1938) with modification (Mann 1945), which requires an observation series y and an accompanying time vector x of length n to detect monotonic changes over time. The Mann–Kendall trend test calculates difference between later-measured data to all earlier-measured data, (y_j—y_i), where j > i are the jth and ith year in the time series, and assigns an integer value of 1, 0, or −1 (positive difference, no difference, and negative difference, respectively) (de Beurs and Henebry 2004). The Mann–Kendall score S is computed as the sum of the integer scores:

$$\begin{equation} S = \mathop \sum \limits_{i = 1}^{n - 1} \mathop \sum \limits_{j = i + 1}^n \left\{ \begin{array}{l@{}l@{}l@{}} {1,\quad \;\;{\rm if}\quad {y_j} - {y_i} > 0}\\ {0,\quad \;\;{\rm if}\quad {y_j} - {y_i} = 0}\\ { - 1,\quad {\rm if}\quad {y_j} - {y_i} < 0} \end{array}. \right. \end{equation} \tag{ 2 }$$

Then the Mann–Kendall test statistic τ is measured by dividing S by the total of n × (n−1)/2 possible pairs of data, where n is the total number of observations for trend direction and strength. Kendall's τ ranges from −1.0 to 1.0, analogous to a correlation coefficient.

We estimated the monotonic rate of change in the time series using the Theil-Sen slope, which computes the slope for all pairs of observations and selects the median value as the robust estimate of the trend's slope (Hirsch et al 1982):

$$\begin{equation} {\beta _1} = {\rm{median}}\left( {\frac{{{y_j} - {y_i}}}{{{x_j} - {x_i}}}} \right). \end{equation} \tag{ 3 }$$

We calculated the area of Theil–Sen slope values associated with a significance level of p < 0.05 for positive and negative trends for each of the seven oblasts of Kyrgyzstan and focused on four rayons (districts). To analyze elevational effects, we divided the area of the focal rayons into five classes: 1400 ≤ × < 1900 m; 1900 ≤ × < 2400 m; 2400 ≤ × < 2900 m; 2900 ≤ × < 3400 m; and × ≥ 3400 m. Note there are no elevations below 1900 m in At-Bashy and Chong-Alay rayons. Note also that Alay and Chong-Alay are adjacent rayons located in the southwest of the country and Naryn and At-Bashy are adjacent rayons in central Kyrgyzstan located to the east of the other pair of focal rayons.

A negative (positive) trend in FDoS indicates earlier (later) onset of snow. A negative (positive) trend in LDoS indicates earlier (later) snowmelt. A negative (positive) trend in DoSS indicates a shorter (longer) snow season.

Table 1. Area in predominant significant (p < 0.05) trends from Terra and Aqua during 2002/3–2015/16 by oblast for snow season metrics: FDoS, LDoS, and DoSS. '–' indicates no prevalent trend.

Oblast	Terra FDoS earlier (km2)	Aqua FDoS earlier (km2)	Terra LDoS earlier (km2)	Aqua LDoS earlier (km2)	Terra DoSS shorter (km2)	Terra DoSS longer (km2)
Batken	526	179	–	–	–	325
Chuy	2079	682	401	375	–	701
Issyk-Kul	–	–	1376	823	884	–
Jalal-Abad	1534	–	759	645	–	–
Naryn	839	–	2227	1242	872	–
Osh	2021	648	–	–	–	701
Talas	972	–	222	213	–	357
Total	7971	1510	4985	3298	1757	2084

To highlight areas of significant change and to attenuate the risk of finding a significant difference where none exists (i.e. a Type I inferential error), we calculated the ratio of area of significant negative trends to the area of the significant positive trends for each administrative unit. We interpreted a factor of >2.0 (or <0.5) in the ratio of significant trends to indicate the predominant direction of change over the study period. We report here only the significant trends showing a predominant direction of change at the level of administrative unit or elevation class. Pixel totals can vary among metrics due to the exclusion of pixels exhibiting no variation in snow cover and thereby generating NaNs (i.e. not a number) in the trend analyses.

Finally, at the national level only, we tracked the trend status (positive at p < 0.05 or negative at p < 0.05 or not significant at p ≥ 0.05) of every pixel for the three metrics at two sequential stages yielding 3² combinations for each of the three metrics as follows:

(1) FDoS = {+ | - | ns} AND DoSS = {+ | - | ns}

(2) LDoS = {+ | - | ns} AND DoSS = {+ | - | ns}

(3) FDoS = {+ | - | ns} AND LDoS = {+ | - | ns}.

Random spatial variation should yield an approximately equal proportion of positive and negative trends. Accordingly, deviations from equal proportions are particularly interesting, especially when the deviations occur in two sequential stages.

We first present the areal extent in each oblast associated with the predominant trend direction and compare the results from the most recent MODIS snow product from Terra (MOD10A2 v006) with the similar MODIS product from Aqua (MYD10A2 v006). Since these satellites have different equatorial daytime overpass times (1030 and 1330, respectively), we expect to see some areal differences. We then focus on the four rayons in the southern zone where we have conducted summer field work in pastures—Alay and Chong-Alay rayons in Osh oblast and At-Bashy and Naryn rayons in Naryn oblast—and compare the areal extent of the predominant trends by elevation class in the Terra product only. (Some tables appear in the supplementary materials available online available at stacks.iop.org/ERL/13/065006/mmedia.)

For visual display only, we identified pixels from Terra with values of Thiel-Sen slope that were significantly different from zero at three significance levels (p < 0.1, p < 0.05, and p < 0.01). The areal analyses are limited to those data significant at p < 0.05.

Earlier snow arrival corresponds to a negative trend in the FDoS, and it appears prevalent across Kyrgyzstan (figure 1(b)), particularly in the oblasts of Chuy (2079 km²), Jalal-Abad (1534 km²) and Osh (2021 km²), for almost 8000 km² in total (table 1, column 1). Only Issyk-Kul oblast does not exhibit a significant predominant trend of earlier snow arrival. Some patches of significant positive trends in FDoS (later snow arrival) are evident in eastern Kyrgyzstan (figure 1(b)), but they are not predominant, resulting in no entry ('–') in table 1.

Earlier snowmelt corresponds to a negative trend in the LDoS, and it appears prevalent across Kyrgyzstan (figure 1(c)), particularly in the oblasts of Naryn (2227 km²) and Issyk-Kul (1376 km²), for almost 5000 km² in total (table 1, column 3). Note, however, that neither Batken nor Osh oblast exhibit snowmelt that is either significantly earlier or significantly later.

A shorter (longer) snow season corresponds to a negative (positive) trend in the DoSS, and the results are mixed across Kyrgyzstan (figure 1(d)). A shorter snow season is apparent in Naryn (872 km²) and Issyk-Kul (884 km²) oblasts, totaling 1757 km² (table 1, column 5). In contrast, a longer snow season appears in parts of four oblasts to the west: Batken (325 km²), Chuy and Osh (each at 701 km²), and Talas (357 km²), totaling 2084 km² (table 1, column 6). Only Jalal-Abad oblast exhibits no predominant change in areal extent of DoSS.

Trend analysis with the Aqua product (MYD10A2) shows predominant trends in fewer oblasts and much smaller areas than in the Terra product (MOD10A2) (table 1, columns 2 and 4). There are no oblasts with predominant positive trends in FDoS, LDoS, or DoSS as well as no predominant negative trends in DoSS (table 1). Earlier snow arrival is evident in just three oblasts in the Aqua product: Chuy (682 km²), Osh (648 km²), and Batken (179 km²) (table 1, column 2). These areal extents are less than the corresponding areas in the Terra product by 67%, 68%, and 66%, respectively. The total area with earlier onset of snow is 1510 km² for Aqua versus 7971 km² for Terra, or a decrease of 81%. Earlier snowmelt (negative trend in LDoS) is evident in five of seven oblasts in the Aqua product, but the proportional decreases in area are less than seen with the FDoS: Chuy (6%), Issyk-Kul (40%), Jalal-Abad (15%), Naryn (44%), Talas (4%), and 34% overall (table 1, column 4).

Based on the Terra data, elevational patterns are consistent across the four focal rayons for LDoS. In both Naryn and At-Bashy, earlier snowmelt occurs at every elevational stratum, with an increase area with increasing elevation most apparent in At-Bashy (table 2, column 1–2). In Alay and Chong-Alay, there are predominant areas of later snowmelt only above 3400 m, but earlier snowmelt below 3400 m (table 2, columns 3–4). The area of earlier snowmelt is greater in the eastern than the western rayons by a factor of 15.

In contrast to the LDoS, the FDoS shows elevational variation across rayons. In Naryn, 167 km² show earlier snow onset and all of it occurs below 3400 m (table S1, column 1). In the neighboring rayon of At-Bashy, there is both earlier snow onset and later snow onset, with the latter appearing between 2900–3400 m (table S1, column 2). Alay exhibits earlier snow onset below 2900 m and above 3400 m (table S1, column 3). Chong-Alay rayon, in contrast, shows earlier snow onset particularly between 2900–3400 m, but not above 3400 m (table S1, column 4).

The DoSS show both longitudinal and elevational patterns of significant change. The two eastern rayons (Naryn and At-Bashy) exhibit areas with shorter snow season duration, with a clear pattern of increasing area with increasing elevation above 2400 m in At-Bashy, but only at the highest elevation class in Naryn (table S2). The two western rayons (Alay and Chong-Alay) show areas of longer snow season duration, with Chong-Alay's duration highest between 2400–2900 m but absent above 3400 m (table S2).

We conducted the two-stage trend analyses only at the national scale and used only the Terra dataset, since it had larger areas of significant change. The varying amount of total pixels in the first stage metrics arises from the exclusion of pixels exhibiting no variation in snow cover, which generated NaNs in the trend analyses.

Two-stage trend analysis for FDoS followed by DoSS shows a substantially larger area with significantly earlier FDoS than significantly later and, of those significantly earlier FDoS pixels, a substantially larger area (19.8%) exhibits significantly longer DoSS than shorter (table 3, row 1). In contrast, there is a substantially larger area (60.1%) showing significantly shorter DoSS associated with pixels with significantly later FDoS (table 3, row 2). No predominant trend in DoSS was evident in those pixels with no significant trends in FDoS (table 3, row 3).

Two-stage trend analysis for LDoS followed by DoSS shows a substantially larger area (3.2%) in significantly earlier LDoS than significantly later and, of those significantly earlier LDoS pixels, a substantially larger area (8.2%) exhibits significantly shorter DoSS than shorter (table S3, row 1). In contrast, there is a substantially larger area (19.9%) showing significantly longer DoSS associated with pixels with significantly later LDoS (table S3, row 2). No predominant trend in DoSS was evident in those pixels with no significant trends in LDoS (table S3, row 3).

Two-stage trend analysis of FDoS followed by LDoS shows a substantially larger area (4.8%) in significantly earlier FDoS than significantly later and, of those significantly earlier FDoS pixels, a substantially larger area (2.3%) exhibits significantly earlier LDoS than later (table S4, row 1). Likewise, there is a substantially larger area showing significantly earlier LDoS associated with pixels with significantly later FDoS (0.9%) or no significant trend in FDoS (91.3%) (table S4, rows 2–3).

Our results, based primarily on the most recent version of Terra MODIS snow cover composites, indicate that snow seasonality has been changing in recent years in each of the seven oblasts of the Kyrgyz Republic: specifically, more territory has been experiencing earlier onset of snow than earlier snowmelt, and roughly equivalent areas have been experiencing longer and shorter duration of snow seasons in the past 14 years (table 1). Significant trends apparent in the Aqua data were less prevalent (table 1). This discrepancy between the Terra and Aqua results may arise from the early afternoon overpass of Aqua, when imaging geometry, cloudiness, and surface temperature may differ. However, the general trend toward earlier snowmelt seen in the Terra product was also evident in the Aqua product.

Zooming into the rayon level and stratifying by broad elevational bands within the focal rayons revealed trend variation in snow cover metrics. We found elevational variation in each snow metric: (1) earlier onset of snow below 2900 m in all four rayons, but divergence above 2900 m, including later onset of snow in At-Bashy between 2900 and 3400 m (table S1); (2) earlier snowmelt below 3400 m in all rayons except Chong-Alay, but later snowmelt above 3400 m in Alay and Chong-Alay (table 2); and (3) occurrence of shorter snow seasons in At-Bashy above 2400 m and longer snow seasons in Chong-Alay below 3400 m (table S2).

At the national level, there are ~8600 km² with significantly earlier snow onset (tables 3, S4) and ~5500 km² with significantly earlier snowmelt (table S3). Regardless of trend status in the first date of snow, there is a substantially larger area exhibiting significantly earlier snowmelt across Kyrgyzstan (table S4). Snow season duration is significantly longer across ~1700 km² where snow onset is significantly earlier (table 3). However, snow season duration is significantly shorter across 980 km² where snow onset is significantly later (table 3). Snow season duration is significantly shorter across 452 km² where snowmelt is significantly earlier (table S3). In contrast, snow season duration is significantly longer in 155 km² where snowmelt is significantly later (table S3). Each of these two-stage trends exhibits strong asymmetry in the area associated with the pair of significant trends, strengthening the interpretation that these trends are not a result of random spatial variation.

Spatial, temporal, and elevational variations in snow seasonality in Central Asia have been detected and quantified in earlier studies at multiple spatial extents. (We summarize the results of the following studies and ours in table S5.) Dietz et al (2013) processed daily MODIS snow cover products between 2000 and 2011 to characterize interannual variation in snow cover across Central Asia with a view to estimating the water content in major regional catchments contributed by snowmelt. They found high spatial and temporal variation in snow cover and no discernable trend in the start, end, or duration of the snow season at this broad scale of analysis.

In a follow-on study, Dietz et al (2014) expanded the temporal scope of the analysis from 1986–2014 by adding coarser spatial resolution Advanced Very High Resolution Radiometer (AVHRR) data to the MODIS data. They divided the snow cover duration into early (01SEP-15JAN) and later (16JAN-31AUG) seasons to detect trends within nine major catchments in Central Asia and by 100 m elevational increments across Central Asia. They found significant positive trends in early season snow cover duration for most catchments, but mixed results for later season snow cover duration, with five of the nine catchments having no significant trend and the remaining four showing significant negative trends. A similar pattern was apparent in their elevational analysis: significant increasing snow cover duration in the early season at most elevations, but more than half the trends in the later season were not significant. Of those that were significant, there was decreasing later season snow duration between 2500 m and 3300 m (Dietz et al 2014). This elevational range is significant: most of the highland pastures on which Kyrgyz agropastoralism depends fall into this band. However, the scale of their analysis was broad, encompassing all of Central Asia.

Zhou et al (2013) used AVHRR and MODIS data to study snow cover trends across the basin of the Amu Darya from 1986–2008. They found statistically significant negative trends in snow cover duration, date of snow cover onset, and date of snowmelt across most of the basin, except trends of earlier snow onset in the Central Pamir, especially at elevations greater than 4000 m (Zhou et al 2013).

Tang et al (2017) found that maps of snow cover days based on cloud-screened MODIS daily snow products exhibited a high mean (> 85%) consistency with in-situ observations of snow cover days. Their analysis focused on the Tien Shan range divided into four regions of which Central Tien Shan corresponds most closely to our study area. Using simple linear regression to identify and characterize trends, Tang et al (2017) found in the mean snow-covered area in each of four seasons in the Central Tien Shan ranging from a minimum of −8% in autumn to maximum of −14% in spring, but no trend was statistically significant at p < 0.05. They also calculated a decrease of −12% in the duration of snow cover in the Central Tien Shan from 2001–2015, but it was also not significant.

Despite substantial spatio-temporal variation within Kyrgyzstan and across the rest of Central Asia, significant trends in snow seasonality exist and these changes have the potential to disrupt herder livelihoods. We took along preliminary maps of the snow season analyses during a field campaign in Osh oblast during July 2017 to sample pastures. The longtime head of the pasture committee for Chong-Alay rayon studied our trend maps and commented that earlier onset of the snow season had indeed become a serious problem for herders in recent years by limiting the use of the fall and winter pastures (Maatkarimov 2017). Note that when we speak of trends, we are talking retrospectively about changes that have already occurred. We are not inferring changes to the local or regional climates using 14 years of data. However, there have been sufficient observations gathered through remote sensing to document significant trends in snow seasonality across large areas in Kyrgyzstan and elsewhere in montane Central Asia.

Recent studies over the Central Asia region show a shift in precipitation from snow to rain. It causes a decrease in snowfall fraction, reducing snow and glacier accumulation during winter (Chen et al 2016). Moreover, changes in snow cover and substantial shrinkage of glaciers induce alterations in the local water cycle, changing runoff and groundwater storage. Less precipitation in the form of snow leads to earlier melting of snow, which can eventually shift peak runoff and river flow toward earlier in the year instead of during the summer when demand for water is highest (Barnett et al 2005, Tang et al 2017). Higher increases of temperature are projected for summer and fall seasons, while lower increases are projected during winter (Xu et al 2017), and a significant decrease in precipitation in spring and summer (Hijioka et al 2014, Yuan-An et al 2013). However, these precipitation projections are highly uncertain (Flato et al 2013, Hijioka et al 2014). The assessment of climate change effects on snow cover is particularly difficult (Tang et al 2017) because those effects strongly vary with geographic context and elevation. Complex terrain generates many local microclimates with different feedbacks making them harder to compare. Sunshine duration, vapor pressure, wind velocity, and their interactions may also enhance spatial and temporal variation in snow cover.

Snow cover affects surface climate, including subsequent vegetation growth (Groisman et al 1994a, Dye and Tucker 2003). Changes in the timing of snow arrival and snowmelt may also have impacts on montane vegetation (Inouye 2000, 2008). Changes in vegetation community composition at higher elevations occur as cold-adapted species decrease in abundance while warm-adapted species increase in a process called thermophilization (Gottfried et al 2012). Furthermore, pasture degradation, including the spread of weedy and unpalatable species, is already a concern in the highland pastures of Kyrgyzstan (Hoppe et al 2016, Eddy et al 2017). A logical next step is to link snow cover seasonality—timing of snow onset, snowmelt, and duration of snow season—with subsequent land surface phenology to detect moisture-induced vegetation stress in highland pastures.

Climate change impacts threaten the prospects for economic development across Central Asia, but especially rural livelihoods that depend on natural resources (Reyer et al 2017). One possible consequence is increased rural to urban migration (Reyer et al 2017). Yet Central Asia in general and Kyrgyzstan in particular are lagging in institutional preparation for climate change adaptation (Ford et al 2015, Lesnikowski et al 2015).

The Working Group II of the IPCC Fifth Assessment report noted that there are manifold knowledge gaps about the impacts of climate change in Central Asia (see table 24-2 in Hijioka et al 2014). This study attempts to help address some of these gaps at a scale finer than most of the literature to date. By using the most recent MODIS snow cover composited product at multiple scales relevant to herder livelihoods in Kyrgyzstan, we have identified areas where snow season timing and duration have already significantly changed in the past 14 years.

Significant shifts toward earlier onset of snow have been identified in nearly 8000 km² in six of seven oblasts and significant shifts toward earlier onset of snowmelt in nearly 5000 km² in five of seven oblasts. In the past 14 years, the duration of the snow season has significantly shortened in two oblasts and significantly lengthened in four oblasts. At finer scales, changes in snow seasonality have varied by elevation and rayon, and the changes detected may impact the montane agropastoralism that forms the basis of the economy in much of rural Kyrgyzstan. The next step is an assessment of how these recent changes in snow seasonality have affected the highland pastures upon which rural livelihoods depend.

We appreciate the assistance and feedback from A A Aidaraliev, K Kelgenbaeva, and P Maatkarimov. We thank the two anonymous reviewers for their useful feedback that helped us to improve the article's clarity. This research was supported, in part, by the NASA Land Cover / Land Use Change Program (NNX15AP81G) How Environmental Change in Central Asian Highlands Impacts High Elevation Communities.