In this study, Central Asia refers to the region comprising the Turkmenistan, Uzbekistan, Tajikistan, Kyrgyzstan and Kazakhstan, as well as northwestern China (Cowan, 2007). Some ice-existing regions in Afghanistan and Pakistan are also included for comparison purpose (Jacob et al., 2012; Yi and Sun, 2014). The ice-covered regions in Central Asia (including those in Afghanistan and Pakistan) are mostly concentrated on the Tianshan Mountains and the Pamir Plateau (Fig. 1). The Tianshan Mountains extend from eastern Uzbekistan to eastern Xinjiang of China with a total length of 2500 km. The total area reaches 390× 10³ km² with 14× 10³ km² covered with ice (glaciers). The Pamir Plateau is actually a mountain knot of three ranges (i.e., Himalaya, Karakoram, and Hindukush). The total area reaches 60× 10³ km² with 30× 10³ km² covered with ice (glaciers). These glaciers feed major rivers over Central Asia (e.g., Tarim, Amu Darya, and Syr Darya), and the variations in glacier mass will thus directly impact the downstream runoff.

Fig. 1

	Figure Option ViewDownloadNew Window
	Fig. 1 Locations of glaciers over Central Asia

As shown in Figure 1, the Randolph Glacier Inventory was used to construct the glacier mask (0.5° × 0.5° grid) over Central Asia (Arendt et al., 2012). The GRACE gravity field spherical harmonic coefficients (SHC’ s, Level 2 (L2) Data Products; Bettadpur, 2003) were used to compute the glacier mass. Monthly GRACE data sets were obtained from the University of Texas’ s Center for Space Research (CSR). Because the atmospheric and oceanic mass changes and tidal effects were presumably removed from GRACE L2 RL05 data, the GRACE gravity field variations primarily reflect the combined effects of ice, snow, waters stored on the surface and in the ground. The variations also contain signal of solid earth effects associated with the ongoing viscoelastic response of the Earth crust and mantle to ice load variability (Chen et al., 2006; Velicogna and Wahr, 2006). We used the GRACE L2 RL05 data covering the period from July 2002 to March 2015 to estimate the glacier mass. It should be particularly pointed out that the GRACE L2 RL05 data were significantly improved over RL04 data in terms of accuracy and also over Level-1B data in terms of data processing techniques (Chambers and Bonin, 2012).

Due to the orbit geometry and the short separation between two satellites (220 km), the time-dependent variable, gravity field coefficients (i.e., C₂₀ value), cannot be reliably determined from the GRACE observations. Thus, we replaced the GRACE-based C₂₀estimate with more reliable value derived from satellite laser range (Cheng and Tapley, 2004). In addition, the SHC’ s (i.e., C₁₀, C₂₀and S₁₁) in degree-1 were cited from the estimates of Swenson et al. (2008) for no provision in GRACE solutions. Also, a decorrelation filter P11M5 (Swenson and Wahr, 2006) was applied to reduce the influence of longitudinal strips, and a 300-km Gaussian filter (Wahr et al., 1998) was further applied to damp the GRACE errors for removing the contamination of SHC’ s in high degrees and orders. The means of all 141 monthly solutions were then removed to obtain time series of gravity field variations (δ C_lm(t), δ S_lm(t)).

To single out of the signal of glacier mass, we used monthly-scale soil moisture estimates at 0.25° × 0.25° resolution from the Global Land Data Assimilation System (GLDAS/Noah) to remove soil moisture-related hydrological effects (Rodell et al., 2004). Besides, a model of glacial isostatic adjustment (GIA) was adopted to remove solid Earth effects (Geruo et al., 2013). The GLDAS and GIA gridded fields were then converted to SHC’ s and the SHC’ s were then truncated into the same degrees and orders with GRACE gravity field. After the Gaussian filter was applied, the SHC’ s of GLDAS and GIA were removed from time series of gravity field variations (δ C_lm(t), δ S_lm(t)).

The processed values of SHC’ s were then used to compute monthly glacier mass over the glacier mask following the method described by Wahr et al. (1998). The missing values of 11 months in the 141 monthly time series of glacier mass were interpolated by least-squares method. Figure 2 shows the patterns in equivalent water height over Central Asia after removing GIA and surface water variations. As shown in Figure 2, the spatial patterns of glacier mass in Central Asia were highly heterogeneous, and an obvious loss of glacier mass occurs in the eastern margin of the Tianshan Mountains.

Fig. 2

	Figure Option ViewDownloadNew Window
	Fig. 2 Spatial patterns in equivalent water height over Central Asia after removing glacial isostatic adjustment (GIA) and surface water variabilities in the Tianshan Mountains (red) and the Pamir Plateau (blue)

As a non-parametric technique and the extension of singular spectral analysis (SSA), multi-channel singular spectral analysis (MSSA) can be used to extract the known and unknown oscillations from noisy spatio-temporal data (Vautard and Ghil, 1989). The extracted oscillations by MSSA can be modulated in amplitude and phase, being consistent with the intrinsically nonlinear dynamics of climate change (Plaut and Vautard, 1994). Given the strong seasonal variability in glacier mass, a 13-month running mean was applied to extract the inter-annual variability. Firstly, the principal component analysis (PCA) was performed on the normalized glacier data set. The first ten principal components (S-PCs) in the PCA analysis were then chosen for the subsequent MSSA study (Plaut and Vautard, 1994; Ghil et al., 2002). According to Vautard et al. (1992), a proper use of a window length M in a time series analysis can successfully isolate the oscillations with certain periodicities and the window length is a critical user-prescribed parameter that defines the spectral resolution of the algorithm. Moreover, the choice of M is a trade-off between quantity of information extracted (a large M) and the degree of statistical confidence in that information (as many repetitions as possible) (Ghil et al., 2002). For a dataset , which is a multivariate time series with L variables (channels) and length T, the starting point of MSSA computational algorithm is to embed each channel into a phase space of dimensionMwith a sliding window of length M (i.e., , t=1, …, T-M+1). Thus, the columns and length trajectory matrix can be obtained. The specific algorithm is as follows:

The block-Toeplitz matrix then can be formed by

Each block is the M× M lag covariance matrix between channel l and channel l’ . As discussed in Plaut and Vautard (1994), the estimate of the entries (l, l’ ) in each block can be calculated as:

Next, diagonalizing the LM× LM matrix gives the eigenvalues and associated eigenvectors (also called space-time empirical orthogonal functions, i.e., ST-EOFs . The eigenvalues are then arranged in decreasing order and help to find the amount of covariance contributed by each mode. Projecting the original seriesXonto the ST-EOFs yields space-time principal components, i.e., ST-PCs.

with k=1, …, LM and n=1, …, T-M+1.

Finally, the original time series are decomposed into reconstructed component (RC) by convolving corresponding ST-EOF with ST-PC (Plaut and Vautard, 1994).

According to Plaut and Vautard (1994), oscillatory components are formed by pairs of ST-EOFs and their associated ST-PCs and can be determined through three fundamental properties of MSSA: (1) two consecutive eigenvalues are nearly equal; (2) two corresponding time sequences described by the ST-EOFs are nearly periodic in the same period and in quadrature; and (3) the associated ST-PCs are in quadrature.

3.1.1 Inter-annual oscillatory pairs

The largest possible window length of 76 was used to extract the inter-annual oscillations. Thus the LM=10× 76=760, and the size of the covariance matrix to be diagonalized is 760× 760. Figure 3 shows the eigenvalue spectra of the first 50 ST-PCs of glacier mass from GRACE data. The error bars in Figure 3 are based on the ad hoc estimate of variance of an eigenvalue given

is the maximum number of degrees of freedom in the time series (Unal and Ghil, 1995). And the confidence level of 95% in Figure 3 is obtained from eigenvalues produced by running MSSA on 100 random monthly data sets (Ghil et al., 2002). As shown in Figure 3, the first eight eigenvalues that pass the test at the 95% confidence level account for 95.6% of the total variance.

Fig. 3

	Figure Option ViewDownloadNew Window
	Fig. 3 The eigenvalue spectra of the first 50 space-time principal components (ST-PCs) of glacier mass from Gravity Recovery and Climate Experiment (GRACE) data

Figure 4 shows the eight corresponding ST-EOFs and ST-PCs, individually or in pairs. The third and fourth ST-EOFs and ST-PCs that represent long-term variations account for 21.8% and 8.8% of the total variance, respectively. Since the ST-EOFs and their associated ST-PCs are well in quadrature, the nearly identical eigenvalues of 1 and 2, 7 and 8 from three potential oscillatory pairs account for 54.5% and 4.3% of the total variance, respectively.

Fig. 4

	Figure Option ViewDownloadNew Window
	Fig. 4 The first eight space-time empirical orthogonal functions (ST-EOFs) and space-time principal components (ST-PCs) of glacier mass

In previous studies, two spurious jumps were found in the atmosphere and ocean de-aliasing level-1b (AOD1B) data, which may cause uncertainties in the estimates of regional glacier mass (Duan et al., 2012). Hence, the L2 GAC data computed from 6-hourly AOD1B data sets were restored to GRACE data to correct the atmospheric effects. Table 1 shows the variabilities of the GRACE glacier mass and the corrected atmospheric effects of the six eigenmodes. The fact that the variance differences between the GRACE glacier mass and the corrected atmospheric effects are within 0.2% seems to suggest that the results of oscillatory pairs was barely affected by the atmospheric effects.

Table 1

Table 1

Table 1 The variabilities of Gravity Recovery and Climate Experiment (GRACE) glacier mass and corrected atmospheric effects of six eigenmodes Eigenmodes Case 1: variance in GRACE glacier mass (%) Case 2: variance in the corrected atmospheric effects (%) 1& 2 54.5 54.3 3 21.8 21.8 4 8.8 8.8 7& 8 4.3 4.4

Table 1 The variabilities of Gravity Recovery and Climate Experiment (GRACE) glacier mass and corrected atmospheric effects of six eigenmodes

3.1.2 Patterns of the inter-annual oscillations

For each channel l, the oscillatory components RC12 and RC78 for the eigenmodes (1, 2) and (7, 8) were obtained by adding the two RCs of each oscillatory pair. Using the maximum entropy method (Burg, 1967), two oscillatory components with periodicities of 6.1-year and 2.3-year were identified from the channellthat has a maximum variance (Fig. 5).

Fig. 5

	Figure Option ViewDownloadNew Window
	Fig. 5 Power spectra from channel l that has a maximum variance in RC12 (6.1-year) and RC78 (2.3-year)

As described in Section 1.3, space-time series of these two oscillatory components(STRCs) can be estimated as:

Where x is the position of grid point over ice-covered regions and t represents time. Figure 6 shows the differences in spatial patterns of the 6.1-year and 2.3-year oscillations between extreme high and extreme low phases, The extreme high and extreme low phases are defined as the maximum and minimum peaks of theSTRCat each grid point. The difference in the 6.1-year oscillation is high in the Pamir Plateau and also in the eastern Tianshan Mountains (Fig. 6a), and, the maximum value (up to 9.6 cm) occurs on the edge of the western Hindukush. However, the difference in the 2.3-year oscillation (Fig. 6b) is much lower in comparison with that in the 6.1-year oscillation and the maximum value (only 2.6 cm) occurs on the eastern Pamir Plateau.

Fig. 6

	Figure Option ViewDownloadNew Window
	Fig. 6 Differences in spatial patterns of the 6.1-year (a) and 2.3-year (b) oscillations between extreme high and extreme low phases

3.1.3 Relationship between inter-annual oscillations and precipitation

A strong correlation between glacier mass balance and precipitation was previously found over Central Asia (Yi and Sun, 2014). To determine the relationship between the inter-annual oscillations in glacier mass and precipitation over Central Asia, we used the monthly precipitation dataset from the Tropical Rainfall Measuring Mission (TRMM) (Huffman et al., 2007). And the precipitation data were processed in the same way as the GRACE data. The TRMM gridded fields were adjusted into the same degree and order as GARCE data. After the same Gaussian filter used on GRACE was applied, the monthly precipitation data were then interpolated to fit into the glacier mask of the Central Asia. When MSSA was performed on the first ten S-PCs of monthly filtered precipitation data, the 6.1-year oscillation was well captured by the first and second eigenmodes, which accounts for 30.4% of the total variance. For the 6.1-year oscillation, the lagged cross-correlation function was computed between MSSA-filtered glacier mass and precipitation on each grid as a function of lead time in precipitation. Figure 7 shows the 6.1-year oscillations for glacier mass and precipitation and corresponding lagged cross-correlation between them at grid point (39.25° N, 67.75° E). The change in the glacier mass was lagged about 13 months behind the precipitation change (Fig. 7a), and the corresponding cross-correlation result shows a significant positive time lagging (Fig. 7b). Figure 8 shows the spatial patterns of the lag time of glacier mass responses to the precipitation over Central Asia. Generally speaking, the glacier mass over Central Asia was lagged behind the precipitation by 9-16 months with an average of 13 months in terms of the 6.1-year oscillation. However, the lagging effect was not detected for the 2.3-year oscillation.

Fig. 7

	Figure Option ViewDownloadNew Window
	Fig. 7 (a) The 6.1-year oscillations for glacier mass and precipitation at grid point (39.25° N, 67.75° E); (b) Corresponding cross-correlation coefficient between the two at grid point (39.25° N, 67.75° E)

Fig. 8

	Figure Option ViewDownloadNew Window
	Fig. 8 Spatial patterns of the lag time of glacier mass responses to the precipitation on each grid (in the 6.1-year oscillation)

We conducted the time series analysis to detect possible trends existed in ST-RC3 and ST-RC4 data (denoted as STRCs34) and ST-RC1, 2, ST-RC3, 4 and ST-RC7, 8 data (denoted asSTRCs123478) of glacier mass over the Tianshan Mountains and the Pamir Plateau (Fig. 9). The results based on the data of STRCs34 only express the general decreasing trends over both the Tianshan Mountains and the Pamir Plateau (i.e., dashed lines in Figs. 9a and b). However, the results based on the data of STRCs12345678 show large-amplitude fluctuations in glacier mass during the period from 2002 to 2015 over both the Tianshan Mountains and the Pamir Plateau (i.e., light-grey lines in Figs. 9a and b). It is clear that inter-annual oscillations are obviously superimposed upon trend in STRCs123478, which can affect the glacier mass trend at short time intervals.

According to published estimates, glacier mass trends over Central Asia inferred from GRACE data differed at different time periods (Matsuo and Heki, 2010; Jacob et al., 2012; Gardner et al., 2013; Yi and Sun, 2014). And this was usually caused by two reasons. Firstly, leakage effect and signal attenuation caused by smoothing and limited spatial resolution made GRACE data problematic to estimate glacier mass trend (Wouters et al., 2008; Chen et al., 2013). Secondly, glacier mass trend from GRACE data might also be affected by the measurement errors in the gravity fields and uncertainty of the background geophysical models used in GRACE data processing (Velicogna and Wahr, 2013). Although the discrepancy of the estimates, the results confirmed the importance of modelling of inter-annual oscillations of glacier mass over Central Asia based on the method.

Fig. 9

	Figure Option ViewDownloadNew Window
	Fig. 9 The trends based on STRCs34 data (dashed line) and STRCs123478 data (light-grey line) in glacier mass over the Tianshan Mountains (a) and the Pamir Plateau (b)