Complex principal component analysis of mass balance changes on the Qinghai – Tibetan Plateau

Climatic time series for Qinghai–Tibetan Plateau locations are rare. Although glacier shrinkage is well described, the relationship between mass balance and climatic variation is less clear. We studied the effect of climate changes on mass balance by analyzing the complex principal components of mass changes during 2003–2015 using Gravity Recovery and Climate Experiment satellite data. Mass change in the eastern Himalayas, Karakoram, Pamirs, and northwestern India was most sensitive to variation in the first principal component, which explained 54 % of the change. Correlation analysis showed that the first principal component is related to the Indian monsoon and the correlation coefficient is 0.83. Mass change on the eastern Qinghai plateau, eastern Himalayas–Qiangtang Plateau–Pamirs area and northwestern India was most sensitive to variation of the second major factor, which explained 16 % of the variation. The second major component is associated with El Niño; the correlation coefficient was 0.30 and this exceeded the 95 % confidence interval of 0.17. Mass change on the western and northwestern Qinghai–Tibetan Plateau was most sensitive to the variation of its third major component, responsible for 6 % of mass balance change. The third component may be associated with climate change from the westerlies and La Niña. The third component and El Niño have similar signals of 6.5 year periods and opposite phases. We conclude that El Niño now has the second largest effect on mass balance change of this region, which differs from the traditional view that the westerlies are the second largest factor.


Introduction
Global sea level rise is causing increasing damage to human coastal developments.Storm tides strike coastal areas more frequently and flood damage is intensifying.The erosion of coasts and coastal lowlands causes beaches to recede.Water in coastal regions becomes polluted and farmlands are under sanitation threats.Seawater absorbs heat and expands, causing an increase in global sea levels (Willis, 2003;Antonov et al., 2005).Rising ocean temperatures accelerate the melting speed of polar ice caps and land glaciers.Part of the meltwater directly (meltwater of polar ice caps) or indirectly (meltwater of glaciers) enters the sea through runoff, contributing to rising sea levels (Nguyen and Herring, 2005;Anny and Frédérique, 2011;Shi et al., 2011;Church et al., 2013).Glacial melting also accelerates the loss of human freshwater resources.
The Qinghai-Tibetan Plateau (QTP) contains numerous glaciers and lakes.The QTP covers an area of 47 000 km 2 and the glaciers are the headwaters of several major Asian rivers.The plateau is notable for its high altitude, large area, and variable climate.For example, the southern and southeastern plateau areas are influenced by the Indian and East Asian monsoon circulations, which bring abundant summer rain.The western portion of the plateau, including the Pamirs, is influenced by the westerlies, which produce dry and rainless areas.The interior of the QTP is less influenced by these circulations and has a continental climate (Yao et al., 2012;Yi and Sun, 2014).Yi and Sun (2014) suggested that the Indian monsoon is much stronger than the westerlies and it can influence precipitation in the Pamirs during winter and summer.Compared with the findings of Yao et al. (2012), Yi and Sun (2014) did not consider the influence of the East Asian monsoon.However, we believe that glacial evolution on the QTP is becoming increasingly complex (Fig. 1), because the El Niño climate pattern has become more frequent and is gradually strengthening.We have evidence indicating that this phenomenon will influence glacier development on the plateau.
Glaciers are sensitive to and provide information on climate change.Their melting process records direct and detailed dynamic change information on local and global climates.Glaciologists and meteorologists reconstruct ancient climates and environmental conditions by analyzing samples taken from the plateau glaciers.These data enable predictions of the response relationships between glaciers and climate change over long time periods and allow forecasting of future climate change (Thompson et al., 2006;Yao and Yu, 2007;Yao et al., 2012).However, for plateaus with sparse human populations, it is difficult to obtain glacier time sequences that have high spatial resolution.
The development of space geodetic technology, especially that of earth observations from space, provides highly precise and continuous observations of glacier mass change and water storage variation in remote regions.These data have significantly improved understanding of the mass balance in polar and Asian alpine regions (Chen et al., 2009(Chen et al., , 2011;;Matsuo andHeki, 2010, 2012;Gardelle et al., 2012Gardelle et al., , 2013;;Jacob et al., 2012;Yao et al., 2012;Gardner et al., 2013;Yi and Sun, 2014;Xiang et al., 2016).In the application of Gravity Recovery and Climate Experiment (GRACE) observation data, their methods are generally similar.After subtracting the signals of the glacial isostatic adjustment (GIA) model and terrestrial water storage model from the GRACE data, residual gravity change can be fully attributed to changes in glaciers.However, the dissimilarity of spatial variation and its causes have not been fully explored.
The change of mass balance in the cryosphere results from interactions between glaciers and the atmosphere at different spatial and temporal scales.Principal component analysis (PCA) is a useful method with which to study the time-varying spatial change of mass balance on the QTP (Fenogliomarc, 2000;Wang and Du, 2000).A great advantage of PCA is that it can describe complicated changes of initial data sets with relatively few variables.However, traditional PCA can detect a standing wave but not advancing waves, because of a lack of corresponding phase information.To overcome this disadvantage, Wallace and Dickinson (2010) developed a complex principal component analysis using the frequency domain (FDPC) method.This performs principal component analysis by calculating the vectors of complex features of a relative spectrum matrix.FDPC is the most common method for studying spatiotemporal transmission characteristics.However, if climate change fluctuates over irregular time intervals and the energy of its principal component is distributed in multi-frequency bands, the spatial change image of every frequency spectrum must be analyzed.In this case, it is inconvenient to use FDPC.Compared with FDPC, complex principal component analysis (CPCA) in the time domain is more appropriate (Horel, 1984).The CPCA method transforms original data and its Hilbert transform into a complex time sequence and conducts principal component analysis by calculating the covariance or complex characteristics vector of the cross-correlation matrix.CPCA is an FDPC method for a full-frequency band.When data sets only have a single frequency, CPCA is equivalent to FDPC.Therefore, CPCA can be used to effectively detect transmitting characteristics, especially when the variance of the principal component is distributed across many frequency bands.
In this study, 153 circa monthly gravity solutions from GRACE Release 05 data were used to reproduce spatial changes of mass balance on the QTP.Then, the main components and corresponding spatial modes and time variation of the mass balance were studied using the CPCA technique.The period of each principal component and its time evolution were examined using wavelet amplitude-period spectrum analysis to explore reasons for the spatial differences of mass balance over the QTP.This analysis helps to clarify the response of mass balance to climate change in the QTB region and may clarify the impacts of glacier melt on water resources, ecology, and environmental disasters.

Data
The variation of earth's gravity field reflects the redistribution of mass inside the earth.Over short time periods (in contrast to geologic time), it can be regarded as mass transfer of the earth's surface and shallow fluids.GRACE was jointly developed by the USA and Germany and it has successfully operated for > 10 years.Its monthly gravity solutions have detected changes of 1 mm geoid fluctuation on a 300 km spatial scale and can monitor gravity field variations caused by changes in hydrology and the cryosphere, glacial isostatic ad-justment, and earthquakes (Ramillien et al., 2006;Chen et al., 2007Chen et al., , 2008;;Velicogna, 2009;Rignot et al., 2011).
The time-varying gravity model used in this paper was the Release-05 (RL05) solutions provided by the Center for Space Research (CSR), University of Texas at Austin.The 153 circa monthly GRACE gravity solutions were taken from January 2003 through September 2015 (∼ 12 solutions are missing).Each solution consists of normalized spherical harmonic (SH) coefficients, to a degree and order of 60.The main enhancements in the new releases are the mean gravity model and corrections of various new background models.Some processing algorithms and parameters were improved, including alignments between the star camera data rate, accelerometer, and K-band system (Bettadpur, 2012).Compared with previous data, the RL05 gravity solutions substantially reduce the stripe noise and are able to monitor 1 mm geoid undulation at the spatial scale of 300 km (Bettadpur et al., 2015;Save et al., 2016).However, at high degrees and orders, GRACE spherical harmonics are contaminated by noise, including longitudinal stripes, and filtering is still needed.In our study, the smoothness priors method (Tarvainen et al., 2002;Zhan et al., 2015) was used to remove noise in the spatial domain.Compared with the Gaussian filter, correlated error filter and the combined filter (Gaussian with 300 km smoothing + correlated error), the smoothness priors method has the advantages of less reduction in signal amplitude at high latitude, preservation of greater detail for short-wavelength components in the result, and less signal distortion at low latitudes.Grid statistics of the filtered field show that the results of the smoothness priors method are most similar to the minimum, maximum, and rms values of the original field (Zhan et al., 2015).

Equivalent water height
According to Wahr et al. (1998), surface mass change can be expressed in the form of surface equivalent water height (EWH) as where ρ e is average density of the earth, a is the equatorial radius, and ρ w is water density.Parameter λ is longitude, θ is colatitude, and P m n (cos θ ) is the nth-degree and mth-order fully normalized Legendre function.Parameter k n is the load Love number.c m n and s m n are normalized SH coefficients.

CPCA
Principal component analysis (PCA) was first formulated in statistics by Pearson (1901).Hotelling (1932) further con-tributed to PCA development.The utility of PCA has been demonstrated in many scientific fields, and it has several alternate names, such as singular value decomposition (SVD) (Golub and Loan, 1996;Mandel, 1982) and empirical orthogonal function (EOF) analysis (Lagerloef and Bernstein, 1988;Kaihatu et al., 1998;Zhang et al., 2004).Eigenvector analysis and characteristic vector analysis are often used in the physical sciences and other fields.PCA (Abdi and Williams, 2010;Helena et al., 2000;Wang and Du, 2000) is a multivariate technique that analyzes a data table in which observations are described by several intercorrelated quantitative dependent variables.Its goal is to extract the important information from the table, represent it as a set of new orthogonal variables called principal components, and display patterns of similarity of the observations and variables as points in maps.Mathematically, PCA depends upon the eigendecomposition of positive semidefinite matrices and SVD of rectangular matrices.
Compared with PCA, the CPCA method (Horel, 1984) introduces phase information and it is a useful method for identifying traveling and standing waves (Pfeffer et al., 2010;Kichikawa et al., 2015).CPCA transforms original data and its Hilbert transform into a complex time sequence and conducts principal component analysis by calculating the covariance or complex characteristics vector of the crosscorrelation matrix.
For the CPCA, a complex observation sequence should first be constructed, which is different from the PCA.For a time-varying observation vector u j (t), its Fourier expansion is u j (t) = ω a j (ω) cos(ωt) + b j (ω) sin(ωt) . (2) In this expansion, j stands for the location of the observation point, t is the observation time, and ω is the Fourier frequency.The constructed complex observation vector U j (t) can be expressed as Here, c j (ω) = a j (ω) According to the definition of c j (ω), Eq. ( 3) can be expanded as The real part of Eq. ( 4) is the original observation sequence and the imaginary part is the Hilbert transform of the real part, which does not change the amplitude of each component of u j (t).However, the phase of each spectral component is advanced by π/2.
The traditional PCA is the principal component analysis of the real observation vector, whereas CPCA analysis is the analysis of the constructed complex vector.After normalization of the complex observation vectors, the average value is subtracted from the complex observation vector of each observation point, and then divided by the standard deviation the complex correlation matrix of the observation point can be expressed as Here r j k represents the multiple correlation coefficients between the j th and kth observation points.CPCA compresses information using the least complex eigenvector e j n of the correlation matrix (Eq.5) and the complex principal component p n (t), because the correlation matrix ( 5) is a Hermitian matrix including n real eigenvalues λ. λ j / n i=1 λ i denotes the contribution percentage of the j th principal component.Observation vector U j (t) can be expressed as the sum of N principal components, where * stands for the complex conjugate, and both complex principal components and complex eigenvectors are orthogonal.The nth complex eigenvector element e j n can be expressed as where e j n indicates the multiple correlation relationship between the j th time sequence and nth principal component.s j n and θ j n are respectively correlative order of magnitude and phase.[• • •] t signifies the mean of times.The time sequence elements of principal components can be expressed as the functional form of amplitude T n and phase n .

Wavelet amplitude-period spectrum analysis
Mass balance on the QTP is under the influence of climate change, and it exhibits unsteady quasiperiodic change.After obtaining the temporal change series of principal components in the area, the time-varying changes of the periods and amplitude (energy) require analysis.We used the wavelet amplitude-period spectrum (Liu, 1999;Liu and Hsu, 2012;Zhan et al., 2003) to analyze its time-frequency information, and chose the Morlet wavelet (Morlet et al., 2012) as the basic wavelet.The wavelet amplitude-period spectrum reflects the time-varying amplitude and period of each periodic term (or standardized periodic term).This means that, in this spectrum, the location of extreme points corresponds to the instant period of a periodic signal (or quasiperiodic term) at that moment, whereas the extreme point value corresponds to the instantaneous amplitude of a certain period signal at that moment.The wavelet amplitude-period spectrum of a time sequence f (t) is defined as where Here, the kernel function ψ(t) is the real part of the Morlet wavelet, δ is a constant, and ω 0 is the frequency parameter, C ψ is a constant, and a and b are scale factors of period and time, respectively.
4 Mass change and its CPCA analysis ) surface mass change field (in units of equivalent water height) was calculated from each GRACE spherical harmonic solutions following Eq.( 1).Then, we filtered each surface mass change field using the smoothness priors method (Tarvainen et al., 2002;Zhan et al., 2015) and interpolated missing data using a spline function at each grid point.GRACE mass rate was then estimated at each grid point using least squares to fit a linear trend, plus annual and semiannual sinusoids to GRACE-derived mass change time series.As fitting results, the amplitude values of annual and semiannual terms are constants, so the calculated trend values contain the contributions from the annual and semiannual trends.The 1 • × 1 • gridded data used here do not improve the resolution of GRACE data.The resolution of the calculated data depends on the degree of the RL05 solutions and the GRACE RL05 solutions are limited by the band-limited nature of GRACE orbit configuration (inclination, altitude, and separation of the twin satellites), with an approximate resolution of around 300 km near the equator (Chen et al., 2017).Relevant information is available from NASA websites (https://grace.jpl.nasa.gov/data/get-data/jpl_global_mascons/).One can also calculate smaller grid data using those solutions but the smaller calculated grid data do not indicate more shortwavelength signals in the results.The accuracy of the calculated data remains 1 mm geoid undulation at around 300 km scale.The accuracy of the calculated grid data depends on the accuracy of the RL05 solution itself (Bettadpur et al., 2015;Save et al., 2016), rather than the size of the grid.the mass variation had no obvious trend.We analyzed mass variation in this area during 2003-2015 using CPCA in order to determine the reasons for mass change.Before the CPCA analysis, data of mass change were filtered, and missing data were interpolated at each grid point.Table 1 shows corresponding eigenvalues of the first five principal components and their contribution percentages to mass change in the area.We used the first three principal components for explanation and description.According to Table 1, the results from CPCA of mass variation over the QTP show that the eigenvalues of the first, second, and third principal components are respectively 82.6516, 25.0562, and 8.6290, and their contribution percentages are 54, 16, and 6 %.Together these explain 76 % of the variation of mass balance in the area.Figure 3b and c depict the temporal evolution of the first principal component and its wavelet amplitude-period spectrum analysis results. Figure 3c shows that the periodic com- ponent that affected the first spatial mode is mainly an annual periodic signal, with relatively stable period and amplitude.
The time-sequence wavelet amplitude-period spectrum results show that the period components of the first spatial mode time sequence are simple and are single annual-period signals featuring steady periods.The result of its wavelet amplitude-period spectrum is the same as the result of the wavelet amplitude-period spectrum of the Indian monsoon indices time sequence (Fig. 3d).
We examined possible relationships between the first principal component and the Indian monsoon indices by calculating their lag correlation coefficient and corresponding 95 % confidence interval (CI) using Monte Carlo hypothesis testing (Table 2).The lag correlation coefficient of the first principal component with the Indian monsoon indices was 0.83, a larger value than the 95 % CI of 0.23.The change of the first principal component lags behind that of the India monsoon indices by 30 days.The values are significantly correlated.From the phase information of mass variation and the correlation analysis, we infer that the first spatial mode in the area is strongly controlled by the Indian monsoon, revealing the influence of the monsoon on rainfall in various areas and its spatial evolution.A branch of the monsoon occurs in the QTP via the AB area and proceeds northward over the Tanggula Mountains with gradually declining energy.It is then blocked by the Qilian Mountains and turns westward, forming a circulation.Another branch proceeds north and enters the Qiangtang Plateau from the central and western parts of the Himalayas.It is obstructed there by the Kunlun and Altun mountains and then progresses westward into the Pamirs.The influence of the Indian monsoon accounts for 54 % of mass balance change on the QTP (Table 1).Based on the time sequence of the spatial mode (Fig. 3b), the Indian monsoon has weakened since 2009 and the monsoon change is the main reason for mass balance change in the area.
Figure 4a shows the second spatial mode and its phase information.This mode is mainly manifested as three mass change zones of southeast-northwest orientation: a positive signal in the southern Karakorum-northwestern India region, two negative signals in the AB area, Qiangtang Plateau (E area), Karakorum, and the southern Qilian Mountains.Red arrows in the figure show phase information of the second spatial mode, which has a disordered direction change.They mainly enter the inland plateau from the southeast and affect its mass balance change.
Figure 4b and c show the temporal evolution of the second principal component in the area and its wavelet amplitude-period spectrum analysis.The wavelet amplitudeperiod spectrum of its time series shows that the periodic component of the second principal component is complicated.It mainly contains a semiannual cycle signal, annual cycle signal, 2-4-year, and 6.5-year cycle signals.The semiannual, annual, and 6.5-year cycle signals have the strongest energy.Energy in the 2-4-year cycle signal is relatively weak, and their energies are all unstable.In comparison with the wavelet amplitude-period spectrum of El Niño evolution in corresponding periods (Fig. 4d), we found that both have 6.5-year and annual cycle signals with consistent phase positions.
We also examined relationships between the second principal component and El Niño by calculating their correlation coefficient and corresponding 95 % CI based on Monte Carlo hypothesis testing (Table 2).Their correlation coefficient was 0.30 compared to the 95 % CI of 0.17.Change of the second principal component lagged that of El Niño by 30 days.This result shows a strong correlation between the two.The spatial phase information and wavelet amplitude-period spectrum data suggest that the second spatial mode in the area is mainly affected by climate change related to the East Asian monsoon and El Niño.Its influence is largely divided into two branches.One branch enters the Qinghai Plateau The Cryosphere, 11, 1487Cryosphere, 11, -1499Cryosphere, 11, , 2017 www.the-cryosphere.net/11/1487/2017/through the Sichuan basin, and the other enters the Qiangtang Plateau through the eastern Himalayas, extends to the northwest of the plateau, reaches the Karakorum, and then turns south.
Figure 5a shows the third spatial mode and its spatial phase distribution information (arrows).The third mode is mainly revealed by the features of two regions, a positive signal in the midwestern area (W of 90 • E) and a negative signal in the region of Linzhi (A area).Mass change in other regions is weakly balanced.The red arrow in Fig. 5a shows phase distribution information of the third spatial mode; its direction shows that the mass change has a clear west-to-east configuration.This indicates that the factors behind the change of this mode came from the western direction.
Figure 5b and c show the time change series of the third spatial mode and its wavelet transform spectrum in the area.The results of the wavelet transform spectrum show that the cycle components of this mode mainly contain semiannual, annual, 2-4-year, and 6.5-year cycle signals.In contrast with the results of the second main component, the energy of the time series of the third spatial mode is mainly concentrated in a 2-6.5-year periodic signal; the annual cycle signal is relatively weak.Except for the 6.5-year signal, the energy of the cycle signals is not stable.The phases of the 6.5-year cycle signals in the second and third main components are in opposition, which suggests that their driving mechanisms are in opposition.
Based on the spatial phase information, we conclude that the third spatial mode is mainly affected by the westerlies and La Niña phenomenon.Its influence can be divided into three branches.One branch moves to the north beyond the Karakorum, enters the Tarim Basin, and then reaches the eastern Qinghai Plateau.Another branch moves east beyond the western Himalayas and enters the Qiangtang Plateau.It then meets the East Asian monsoon around 90 • E and is obstructed.The third branch goes southward along the Himalayas and influences northern India.The westerlies are weak in the south and strong in the north, so a clear northeast-southwest boundary of force range (blue line in Fig. 5a) is formed in the inland part of the QTP.

Mass change in inland QTP
In the inland part of the QTP, there are three regions of mass increase: the Qiangtang Plateau (E area), central and eastern parts of the Kunlun Mountains (F area), and Qinghai Plateau (G area).Their respective annual increases were estimated to be 4.5, 5.5, and 3.5 GT and these were much smaller than the 30 GT of Yi and Sun (2014).Many scholars have conducted related research in an attempt to explain the reason behind mass balance change in the region.
Mass balance of the inner Tibetan Plateau (ITP) derived from GRACE data showed a positive rate that was attributable to glacier mass gain, whereas those same glaciers evaluated in other field-based studies showed an overall mass loss.Jacob et al. (2012) deduced glacier mass balance using GRACE data and found a mass increase rate of 7 Gt yr −1 in areas E and F. Yao et al. (2012) observed more than 20 glaciers in the QTP area and concluded that they are shrinking dramatically.Their results indicate that the Himalayas have shown the most extreme glacial shrinkage based on reductions of glacier length and area.The shrinkage is most pronounced in the southeastern QTP, where glacier length decreased at a rate of 48.2 m yr −1 and the area declined at a rate of 0.57 % yr −1 during the 1970s-2000s.The rate of glacial shrinkage decreased from the southeastern QTP to the interior.Zhang et al. (2013) studied 53 % of the total lake area on the plateau using ICESAT satellite data and found a mass increase rate of 4.95 Gt yr −1 .They suggested that the increased mass measured by GRACE was due to increased water mass in lakes.If this rate holds true for all lakes, the total mass variance rate, using the area ratio, is +8.06 Gt yr −1 .However, glacier melting into lakes, by itself, should not increase the overall mass and may decrease the mass because a portion of the meltwater would be lost through evaporation or through discharge into rivers that leave the Tibetan Plateau.Yi and Sun (2014) suggested a relatively large mass rate change in this area, and explained that this change was caused by glacier changes, lake water levels, geologic structural processes, and frozen soil.They stated that, based on model calculation, the change of inland water storage was −3.3 Gt yr −1 .The change of negative balance of weakening glacier mass has been confirmed (Bolch et al., 2010(Bolch et al., , 2012;;Yao et al., 2012).According to calculations of Zhang et al. (2013), the increase of lake water is 8.1 Gt yr −1 , and the effect of tectonic movement (using simple Bouguer correction) is 0-13 Gt yr −1 .The effect of other factors is nearly zero.However, we still lack sufficient observation data of mass balance states in the interior part of the earth in the study region.Thus, the exact Bouguer equilibrium correction requires more data for confirmation.
The effect of soil freezing on mass change in the inland plateau is weak, because the terrain there is flatter than at the plateau edge.The inland area contains numerous lakes and wetlands, which are conducive to fluid convergence.Moreover, when water melts and is lost from frozen soil, soil porosity increases, which captures more water during rainy periods.
Our results indicate that rainfall is the main reason for the mass increase in the study region.There is strong evidence that precipitation over the inland QTP during the past several decades has greatly increased (Yao et al., 2012; Global Precipitation Climatology Project or GPCP, www.esrl.noaa.gov/psd/data/gridded/data.gpcp.html).Through the influence of El Niño, moist air moves westward to the inland plateau through the eastern Himalayas and Qinghai, and this brings rainfall to the inland areas and causes rainfall accumulation in plateau lakes and wetland areas.The inland Plateau, especially the western part of Qiangtang plateau and Kunlun Mountains area, is also influenced by the westerlies and the La Niña phenomenon (Fig. 5a), which further creates the meteorological conditions for rain and snow.Increased temperatures (Qin et al., 2009) accelerate glacial melting in this area.This glacier meltwater enters lakes through runoff.It also explains why on-site observation data of glaciers indicate slight shrinkage, and GRACE observations indicate the reasons for mass increase.

Mass change of glaciers in Himalayan region
The trend of mass balance change from GRACE data shows that the most negative signal is along the Himalayas and northwestern India.The mass reduction rate of glaciers in the entire Himalayan mountain region is 14 Gt yr −1 , and the mass loss of glaciers in the eastern Himalayas was the most dramatic, with the rate of −4.6 Gt yr −1 in the A area and −4.1 Gt yr −1 in the B area.The mass reduction rate in northwestern India (H area) was −13.6 Gt yr −1 , whereas Rodell et al. (2009) and Yi and Sun (2014) estimated larger values of −17.7 and −20.2 Gt yr −1 , respectively.The reason for this difference is that Rodell et al. (2009) used the data of the earlier RL04 version.Yi and Sun (2014) stated that the RL04 solutions overestimate the glacier melt rate in the Himalayas by as much as 17 %.The difference between our results and those of Yi and Sun (2014) stems from their use of the mascon inverse method in a concise form.Moreover, the filtering method used may attenuate the signal.Yao et al. (2012) studied glacial change over the past three decades and found that the Himalayas had the most extreme glacial shrinkage based on reductions of both glacier length and area; the shrinkage was greatest in the southeastern QTP (A area), where the length decreased at a rate of 48.2 m yr −1 and the area was reduced at a rate of 0.57 % yr −1 .Most negative mass balances occurred along the Himalayas and ranged from −1100 to −760 mm yr −1 .This trend of mass change along the Himalayas is consistent with our results.They attributed this change to the weakened Indian monsoon towards the plateau interior.Thakuri et al. (2014) studied glacier changes on the southern slope of Mt.Everest from 1962 to 2011 (400 km 2 ) using optical satellite imagery.They concluded that the observed glacier shrinkage, upward altitude shift of the snowline, and the negative mass balance (Nuimura et al., 2012) are not only due to warming temperatures, but are also the result of weakened Asian monsoons.Bolch et al. (2011) examined the mass change of glaciers on Mt.Everest using stereo Corona spy imagery (1962 and 1970), aerial images, and high-resolution satellite data (Cartosat-1).They found that glaciers south of Mt.Everest continuously lost mass from 1970 to 2007, but at an increased rate in recent years.Wagnon et al. (2013) ar-rived at the same conclusion and noted that glacier shrinkage south of Mt.Everest was less than shrinkage of other glaciers in the western and eastern Himalayas and southern and eastern Tibetan Plateau.Salerno et al. (2015) analyzed the precipitation time series during 1994-2013 reconstructed from seven stations at elevations between 2660 and 5600 m.They found that precipitation has decreased by 47 % during the monsoon period and snowfall has decreased by 10 % in the last 20 years.Salerno et al. (2016) extended this analysis back to the early 1960s and for all regions used, as proxy of the precipitation trend, the surface area variation of glacial lakes.These authors determined that an increase in precipitation occurred until the mid-1990s followed by a decrease until recent years in all Mt.Everest regions.
Studies using different types of data have produced similar results, i.e. negative mass balances and a weakened Indian monsoon along the Himalayas.Our results support these conclusions.The results of CPCA analysis indicate that mass change on the Himalayas and its southern portion are associated with the Indian monsoon climate, and the intensity of this monsoon is weakening.This result is also consistent with the conclusions of Wu (2005).A weakened Indian monsoon brings less humid air to the study region resulting in decreased annual rainfall (Thakuri et al., 2014;Salerno et al., 2015Salerno et al., , 2016)).The GPCP rainfall data confirms this conclusion.The eastern Himalayas are also affected by El Niño (Fig. 4a) and East Asian monsoons, and no evidence supports the role of westerlies (Fig. 5a) in driving local climate and glacier changes.Glaciers in this area are of a marine type, with masses having large inputs and outputs and strongly affected by changes of the marine climate.The weakened Indian monsoon, strengthening El Niño, and westerlies, combined with the huge topographic landform, exert climatic controls on the distribution of existing glaciers along all Himalayan regions and reduce precipitation there.

Effect of circulation in QTP area
Archer and Fowler ( 2004) indicated that the western Hindu Kush and Karakoram are largely exposed to the arrival of westerly midlatitude perturbations bringing precipitation during winter and early spring, whereas the eastern Himalayas is dominated by summer monsoon precipitation (Syed et al., 2006;Yadav et al., 2012).Their results are similar to those of this study.The results of CPCA indicate that the eastern Himalayas is under the influence of weakened Indian monsoon and El Niño, while the Hindu Kush and Karakoram area is under the influence of a weakened Indian monsoon, westerlies, and La Niña.Thompson et al. (2000) examined the variability of the South Asian monsoon by analyzing ice core records of the Dasuopu glacier on the QTP.They found evidence of drought conditions and a weak monsoon from 1780 to 1810.Interestingly, according to historical records, at least 600 000 people died in 1972 in just one region of northern India from an epic drought associated with this event.The onset of this event in the Dasuopu cores is concurrent with a very strong El Niño-Southern Oscillation (ENSO) from 1790 to 1793, which was followed by a moderate ENSO event of 1794-1797.These data suggest an association between ENSO and a weakened Asian monsoon.
Arctic amplification may impact midlatitude weather patterns and extremes (Francis and Vavrus, 2012;Screen and Simmonds, 2013), and midlatitude westerlies may increase climate variation and glacier variability in monsoon affected areas of High Asia (Thomas et al., 2014).On large spatial scales, climate change over the QTP may also be connected with hemispheric or global atmospheric circulations including the North Atlantic Oscillation (NAO) and ENSO (Wang et al., 2003).ENSO may influence climate over the southern QTP through linkage with the Indian monsoon (Xu et al., 2010(Xu et al., , 2011)).The NAO is associated with climate fluctuations over the northern QTP through modulation of the westerlies (Wang et al., 2003;Xu et al., 2010), which is similar to climate change from the westerlies and La Niña in the third principal component.Yao et al. (2012) studied glacial changes over the last three decades in the QTP and found that glacier recession in the Himalayas was the most dramatic, followed by recession in the inland plateau.Glaciers in the Pamirs had weak balance changes, and some of the glaciers in the plateau of the eastern Pamirs continue to expand.Yao et al. (2012) concluded that the main reason for these changes was the variation of climates with different circulations, which includes effects of the weakened Indian monsoon in the Himalayas, rainfall decreases, effects of the strengthening westerlies in the Pamirs and its eastern portion, and rainfall increases.In the inland plateau, the influences of these two circulations are limited.The two atmospheric circulation patterns, combined with the huge topographic landform, exert climatic controls on the distribution of existing glaciers.The East Asian monsoon only affects glaciers on the eastern margin, such as the Mingya Gongga and those in the eastern Qilian Mountains.The interior of the QTP is dominated by continental climatic conditions, and the sparse glacier distribution and higher ELAs in the continental-climate-dominated interior are consequences of a limited water vapor source from both air masses.They divided glaciers of the Tibetan Plateau into seven regions and categorized them into three climatic transects: (1) southwest-northeast oriented (central Himalayas-Qiangtang Plateau-eastern Qinghai Plateau region), with the weakened Indian monsoon influence northward; (2) southeast-northwest oriented (eastern Himalayas-Qiangtang Plateau-Pamirs region), with the weakened Indian monsoon toward the interior and strengthening westerlies toward the northwest; (3) along the Himalayas, with stronger monsoon influence in the east and weaker monsoon influence in the west.
From the results of the CPCA, the first spatial mode shows that the mass balance of the Himalayas-Pamirsnorthwestern India region (transect 3) was the most sensitive to climate change associated with the Indian monsoon, whereas there was little impact of that change on mass balance of the inland plateau.The third spatial mode shows that mass balance of the northwest plateau, including the Kunlun Mountains, is also affected by climate change from the westerlies and La Niña.Another difference between the results of Yao et al. (2012) and those of this study is that climate change from El Niño rather than the weakened Indian monsoon toward the interior affected mass balance along transect 2. We found that the time evolution of the second principal component and El Niño index had a stronger time-frequency correlation.
Yi and Sun (2014) used harmonic analysis of a time series of mass changes in the study region and found a 5-year periodic signal in the Pamirs and Karakorum regions.They analyzed the correlation between mass change, precipitation, El Niño-Southern Oscillation (ENSO) and the Arctic Oscillation (AO), and found that the 5-year undulating signal of mass change is controlled by both the ENSO and AO.The DKMD is located on the eastern Qiangtang Plateau (the center of transect 2), where area mass balance change is the most sensitive to El Niño in our results.Yao et al. (2012) considered the effect of the Indian monsoon and westerlies but did not consider El Niño, which was the second major component (16 %) in the study region.Yi and Sun (2014) noted that the 5-year periodic signal in the Pamirs region is related to ENSO, but ignored the effect of La Niña because they did not distinguish the phase information.According to the CPCA, we believe that the mass change in the QTP area is mainly controlled by the Indian monsoon and westerlies but the influence of El Niño and La Niña on the inland areas of the plateau and the Karakorum area is also important.The Indian monsoon mainly affects mass balance change on southern and southwestern QTP, whereas El Niño mainly modifies that change over the eastern Himalayas, Qiangtang Plateau, Pamirs, and eastern Qinghai Plateau area.Mass balance over the western and north-western QTP is mainly affected by the westerlies and La Niña.

Conclusions
During 2003-2015, mass changes on the Tibetan Plateau and surroundings varied systematically from region to region.Specifically, the Himalayas region had the greatest negative mass balance with a mass decrease at a rate of −14 Gt yr −1 .The continental interior of the plateau had a positive signal with a mass increase at a rate of 13.5 Gt yr −1 .The Pamirs had a weak negative mass balance.The main cause of the systematic mass change was variation in rainfall which mainly resulted from changes in four different atmospheric circulation patterns over the QTP and its surroundings.These were the weakening Indian monsoon, strengthened westerlies, El Niño, and La Niña.Their contributions explained approximately 76 % of mass changes on the QTP.
Change of the Indian monsoon had the most important effect on mass balance variation over the QTP.The lag correlation coefficient of the first principal component with the Indian monsoon indices was 0.83 and much larger than the 95 % CI of 0.23, and the change of the first principal component lags that of the India monsoon indices by 30 days.Mass balance variation over the eastern Himalayas, Karakoram, Pamirs and northwestern India was most sensitive to changes of the Indian monsoon, and was responsible for 54 % of that change.The weakened Indian monsoon, combined with the huge topographic landform, exerted climatic control on the distribution of existing glaciers in these regions and caused less precipitation there.
Because El Niño has been strengthened, it has become the second most important major effect on mass balance change of QTP and is responsible for 16 % of the change.Their lag correlation coefficient is 0.30 compared to a 95 % CI of 0.17 and change of the second principal component lags that of El Niño by 30 days.Mass balance over the eastern Himalayas, Qiangtang Plateau, Pamirs, and eastern Qinghai Plateau areas were the most sensitive to El Niño variation.Further research will increase our understanding of the physical mechanisms linking El Niño and mass balance.
The third principal component was climate change of the westerlies and La Niña.Mass balance on the western and northwestern QTP was the most sensitive to climate change from the westerlies and La Niña, which represented 6 % of the mass balance change.The strengthening westerlies and La Niña climate phenomenon created meteorological conditions conducive for rain and snow to those regions, and our results do not support the role of westerlies in driving glacier changes across the southeastern QTP.

Figure 1 .
Figure 1.Distribution of glaciers (white dots) and atmospheric circulation in and around the Tibetan Plateau.
Figure 2 shows the trend of mass balance on the QTP during 2003-2015.The QTP mass balance has two major change characteristics, namely, a large negative signal with mass decrease around the southern edge of the plateau (Himalayas and its southern region) and a positive signal with mass increase over inland areas of the plateau.In the Pamirs region,
Figure 3a shows the first spatial mode and its spatial phase distribution (arrows) from the CPCA analysis of the mass balance change in the area.The first spatial mode shows change characteristics of three areas: two negative signals from the eastern Himalayas to the Hengduan Mountains (AB area) and the Pamirs to the Karakorum (D area), and a positive signal in northwestern India (H area).The direction of the arrows indicates the sequence of mass change and arrow size indicates the change rate of mass.The phase information demonstrates that the first spatial mode mainly reflects the south-to-north character of mass change.

Figure 3 .
Figure 3. First spatial mode and phase (red arrows) (a), temporal patterns of first principal component (b), and its wavelet amplitudeperiod spectrum (c) of mass balance change, as well as wavelet amplitude-period spectrum of Indian monsoon indices in the period 2003-2009 (d).

Figure 4 .
Figure 4. Second spatial mode and phase (red arrows) (a), temporal patterns of second principal component (b), and its wavelet amplitude-period spectrum (c) of mass balance change, as well as wavelet amplitude-period spectrum of El Niño during the period 2003-2015 (d).

Figure 5 .
Figure 5. Third spatial mode and phase (red arrows) (a), temporal patterns of third principal component (b), and its wavelet amplitudeperiod spectrum (c) of mass balance change.
Ke et al. (2017) examined area and thickness changes of glaciers in the Dongkemadi (DKMD) region of the central QTP using Landsat images from 1976 to 2013 and satellite altimetry data from 2003 to 2008.They analyzed relationships between glacier variation and local and macroscale climate factors based on remote sensing and reanalysis data.Their results indicate that glacier change in the DKMD region was dominated by variation of mean annual temperature, and was influenced by the state of the NAO over the past 38 years.The mechanism linking climate variability over the central QTP and state of the NAO is most likely via changes in strength of the westerlies and Siberian High.In addition, ENSO may have been associated with extreme weather (snowstorms) in October 1986 and 2000 which might have led to substantial glacier expansion in the following years.

Table 1 .
Eigenvalues and contribution percentages to mass change in CPCA analysis of Qinghai-Tibetan Plateau.

Table 2 .
Correlation analysis based on Monte Carlo hypothesis testing.Time lag (month) First principal component Second principal component 95 % confidence level