Southeast Tibet: implication for substantial meltwater contribution to the Brahmaputra

. High Asia glaciers were observed to be retreating the fastest in the southeastern Tibet Plateau (SETP), where vast amounts of glacier and snow (GS) feed the streamflow of the Brahmaputra, a transboundary river linking the world's two most 15 populous countries China and India. However, the low temporal resolutions in previous observations of GS mass balance obscured the seasonal accumulation/ablation variations, and their modelling estimates were divergent. Here we use monthly satellite gravimetry observations from August 2002 to June 2017 to estimate GS mass variation in the SETP. We find that the “spring-accumulation type” glaciers and snow in the SETP reach their maximum in May. This is in stark contrast to seasonal variations in terrestrial water storage, which is controlled by summer precipitation and reaches the maximum in August. These 20 two seasonal variations are mutually orthogonal and can be easily separated in time-variable gravity observations. Our results show a long-term trend of -6.5 ± 0.8 Gt yr -1 (or 0.67 ± 0.08 w.e. m/yr) and annual mass decreases ranging from -43.4 Gt to -73.1 Gt with an average of -57.6 Gt in the SETP between August 2002 and June 2017. The contribution of summer meltwater to the Brahmaputra streamflow is estimated to be 51 ± 9 Gt. This result could help to resolve previous divergent modelling estimates and underlines the importance of meltwater to the Brahmaputra streamflow. The high sensitivity between GS melting 25 and temperature on both annual and monthly scales suggests that the Brahmaputra will suffer from not only changes in total annual discharge, but also an earlier runoff peak due to the ongoing global warming.


Introduction
The Tibetan Plateau, considered as the Asian water tower, is the source of several major river systems. Their upper streams are fed by rainfall, base flow and widespread glaciers and snow (GS) melt (Barnett et al., 2005;Immerzeel et al., 2010;30 Jansson et al., 2003;Lutz et al., 2014). The GS melt is susceptible to climate change, whereas its sustainable supply is critical to the local freshwater security, flood prevention and control, and hydroelectric development (Bolch et al., 2012;Kaser et al., 2010;Yao et al., 2012). The southeastern Tibet Plateau (SETP), including the Nyenchen Tonglha Mountains (NTM) and eastern Himalayas, holds 10,439 glaciers with a total area of 9,679 km 2 (RGI Consortium, 2017) and widespread seasonal snow coverage of up to 100,000 km 2 . These maritime glaciers are characterized by low equilibrium-line altitudes with large 35 topographic gradients (Yao et al., 2012) and the most severe mass loss in High Mountain Asia (HMA) (Brun et al., 2017;Kääb et al., 2015). The GS melt serves as an essential water supplier for the Brahmaputra river system (e.g., Immerzeel et al., 2010;Lutz et al., 2014), which runs through three densely populated countries, China, India and Bangladesh (Figure 1). The revealed vulnerability of glaciers in the Brahmaputra Basin to global warming and emerging controversies over water allocation (e.g., dam building (Tanck and Fazani, 2010)) are increasingly attracting scientific and public concerns. 40 Due to the lack of observational data, most of the previous estimates on the contribution of seasonal meltwater to the upstream flow of the Brahmaputra River were based on modelling approaches that were only calibrated by employing streamflow data. Resultingly, the previous estimates disagree widely from 19% to 35% (Table 1) due to different forcing data and approaches without direct constraints on GS mass balance (Bookhagen and Burbank, 2010;Chen et al., 2017;Huss et al., 2017;Immerzeel et al., 2010;Lutz et al., 2014;Zhang et al., 2013). The amount of meltwater could be even more divergent. 45 For example, Huss et al. (2017) estimated that the amount of annual GS melt to the Brahmaputra River was 138 w.e. (water equivalent) km 3 yr -1 , which is however triple the estimate of 43 w.e. km 3 yr -1 by Lutz et al. (2014). Although these two studies covered different ranges, the glacierized zone in the basin was both well enclosed so the estimates should not be so different.
Such huge discrepancies in previous estimates make it imperative to incorporate the calibration from GS mass balance observations into future modelling experiments. Actually, the concept of assimilating more GS observations has begun to be 50 implemented in the state-of-the-art models (Wijngaard et al., 2017;Biemans et al., 2019), but their glacier results suffered from coarse temporal resolution (two observations over 5 years) and the snow mass changes were partially constrained by area changes.
Spaceborne sensors can be helpful in this desolate mountain region. Remote sensing techniques for region-wide GS mass balance measurements can be divided into three categories: laser altimetry (e.g., Ice, Cloud and land Elevation Satellite 55 (ICESat) (Kääb et al., 2012)), multi-temporal digital elevation models (e.g., SPOT (Gardelle et al., 2013), ASTER (Brun et al., 2017)), and space gravimetry (Gravity Recovery and Climate Experiment (GRACE) (Matsuo and Heki, 2010;Yi and Sun, 2014)). The first two geodetic approaches require the average ice density to convert volume changes into mass changes. The ICESat observation suffers from short operation period (2003)(2004)(2005)(2006)(2007)(2008)(2009) and sparse spatial sampling, both of which can be overcome by the stereo-imagery approach, which is becoming popular for the whole HMA study recently (Brun et al., 2017;60 Dehecq et al., 2018). Brun et al. (2017) provided an estimate of the detailed glacier mass balance trends over HMA between 2000 and 2016 and highlighted the regional dissimilarity. Despite recent improvements in spatial resolution in HMA glacier mass change studies, there has been little advance in their temporal resolution.
Observations at a monthly temporal resolution are necessary to separately quantify summer and winter mass balances, two processes dominating the annual glacier mass balance (Cogley et al., 2011), and thus crucial for the calibration and 65 validation of glaciological models. The amplitude of seasonal variation of the glaciers in the SETP is up to ~3 m w.e. (Wang et al., 2017), far exceeding their net annual change of ~ 0.6 m w.e. (Brun et al., 2017). Hence, the long-term trend of GS mass changes only reflects a small net imbalance of their ablation and accumulation. High time-resolution monthly observations by GRACE since its launch in 2002 (Tapley et al., 2004) are promising in identifying these two processes. Up to now, the application of GRACE in HMA glaciers has been focusing on their secular changes with little attention to the seasonal 70 variations (Gardner et al., 2013;Matsuo and Heki, 2010;Yi and Sun, 2014). This is mostly due to the poor spatial resolution of GRACE (> 300 km) and the dominance of terrestrial hydrological signals in the seasonal gravity signals, which is difficult to eliminate from glacial signals. The latter is particularly severe in the SETP with intense monsoon precipitations. The GS and hydrological mass changes (mainly including mass changes in rivers, soil moisture and groundwater) dominate the seasonal gravity signals in the SETP observed by GRACE. Despite the general difficulty in separating them in the spatial 75 domain, we find it possible to separate the two signals in the time domain, owing to their contrasting seasonal behaviours.
Precipitation in the SETP is controlled by various atmospheric circulation systems in different seasons, with westerly winds and Bay of Bengal vortex in winter/spring and Indian monsoon in summer (Wu et al., 2011;Yang et al., 2013;Yao et al., 2012). The former two systems were found to drive the spring precipitation in the SETP along the Brahmaputra River, thus forming a 'spring-accumulation' type of glaciers . The Indian monsoon prevails from June to September 80 and brings intense precipitation on the southern side of the Himalayas, where terrestrial water storage shows tremendous seasonal changes and peaks in late summer. Therefore, according to the climate stations near NTM, we can observe bimodal precipitation variations throughout the year .
In this work, we will first introduce the precipitation characteristics in this region by both meteorological stations and global precipitation products. We will then use the empirical orthogonal function (EOF) analysis to decompose hydrological 85 and GS signals in our study region, which does not exactly coincide with the range of glacierized zone in the Brahmaputra Basin. Our study region covers only 83% of the basin glaciers (the 17% undetected ones are in the western part) and 15% of non-Brahmaputra glaciers. We will scale our results by a ratio of 1 × #.%& #.%' = 1.02 to get the total meltwater in the Brahmaputra, assuming that our observations can represent the basin-wide average. The hydrological and GS signals are further compared to the results of other datasets to validate their physical meanings. Such high time-resolution observations also allow us to 90 compare GS mass variations with temperature records during the ablation season, and to study the sensitivity of GS mass change in response to the temperature change. Finally, we will compare our results to previous estimates at monthly, annual, and interannual scales.

GRACE data and preprocessing 95
We adopt the monthly GRACE spherical harmonics Release 06 products from August 2002 to June 2017. The three datasets are solved respectively by three organizations: Center for Space Research (CSR) at the University of Texas, ftp://podaac.jpl.nasa.gov/allData/grace/L2/. The degree 1 terms, which are absent in original GRACE releases, have been added based on the technique proposed by Swenson et al. (2008). The C20 terms have been replaced by those from satellite 100 laser ranging (Cheng et al., 2011), which are considered to be more reliable. A widely used Glacial Isostatic Adjustment (GIA) model by A et al. (2013) is adopted to correct the GIA effect caused by historical polar ice sheet changes.
Two different filtering strategies, a combination of P4M6 decorrelation (Swenson and Wahr, 2006) and 300km Gaussian filter (hereafter short for G300+P4M6) and a DDK4 filter (Kusche et al., 2009), are applied separately. Therefore, there are six combinations and their average values (with uniform weights) are used in the following figures. 105

GRACE error estimation
We adopt different uncertainty estimation strategies for the seasonal variation and the trend due to their intrinsically different error sources. The error of seasonal variation consists of the standard deviations among these six datasets (i.e., errors from the data solution and smoothing methods) and the leakage error, while that of the long-term trend also includes other potentially uncorrected signals. The leakage error is determined by how effectively the hydrological and glacial signals are 110 separated by the EOF technique. Based on the modelled and recovered glacier mass changes, their residuals are estimated to have a seasonal variation of up to 11% of the modelled glacier mass change (refer to section 3.2 in the supporting materials), which is used to calculate the seasonal leakage error.
For the long-term trend error, the three different solutions and two smoothing techniques have a total effect of 0.44 Gt/yr.
There are potential errors from other signal sources, like glacial isostatic adjustment (GIA), little ice age (LIA) and weather 115 denudation. The GIA effect which originates from the polar regions has been corrected by A's GIA model (A et al., 2013), although its influence on the trend is as small as 0.02 Gt/yr. The main reason is that the spatial pattern of GIA is quite smooth, so it mainly influences the first mode and rarely leaks into the second one. This feature is also applicable for other signal sources: unless they exactly locate in the glacierized area, their influence will be reduced by the EOF decomposition. In the southern and southeastern Tibetan Plateau (over 500,000 km 2 ), the effects of LIA and denudation are estimated to be -1 ± 1 120 Gt/yr (Jacob et al., 2012) and 1.6 Gt/yr (assuming the sediment has a density of 2 Gt/km 3 ) (Sun et al., 2009), respectively. Our glacierized zone and surroundings have an area of about 100,000 km 2 , accounting for one-fifth of the whole region, so we suppose their contribution to the GS mass estimate is also proportionally 1/5. However, as we explain above, we could not precisely quantify their contribution without knowing their spatial distribution, and they are more likely to be absorbed by the first mode, so we only include their contribution in the error estimation rather than correcting them in the trend. Table 2  125 summarizes the sum of GRACE error estimates in the secular trend.

ICESat altimetry
Version 34 of the ICESat Global Land Surface Altimetry Data is used to derive glacier height changes. The data span is from 2003 to 2009, with two or three observation campaigns per year ( Figure S1). The processing of ICESat data includes the 130 following steps. (1) Orthometric heights are obtained from original elevation data based on the Earth's gravity model 2008.
(2) Footprints on glaciers are identified based on RGI 6.0 glacier outlines. (3) For each ICESat footprint, SRTM (Farr et al., 2007) elevations and slopes are extracted by bilinear interpolation of the DEM grid cells. Glacier height variation is defined as the elevation differences between the footprints and the SRTM data. (4) We exclude footprints over SRTM voids, footprints with slopes higher than 30°, and footprints with height change larger than 100 m (which are attributed to biases caused by 135 cloud cover during the ICESat acquisition). (5) We also discard the calibration campaign L1AB (March 2003) and the incomplete campaign L2F (October 2009). (6) Glacier height variations are averaged and interpolated along the altitude to alleviate the uneven sampling problem in space, and an uncertainty of 0.06 m/yr (Kääb et al., 2012) is chosen to account for the uneven sampling bias in time. The steps have been used in previous work (Wang et al., 2017) and have also been described in earlier studies (Gardner et al., 2013;Kääb et al., 2012). The footprint information is given in Figure S2. 140 ICESat has shown good ability to solve snow variation in flat regions , but applying the same technique in mountainous areas with high terrain heterogeneity is cumbersome. Therefore, here ICESat is only used to estimate changes in glacier mass. Although our GRACE estimate includes both glaciers and snow, the estimates by GRACE and ICESat are comparable in the late ablation season (i.e., the October/November campaign of ICESat), when the contribution of seasonal snow meltwater is negligible. To convert the glacier thickness changes into mass changes, two parameters are required, i.e., 145 glacier density and total glacier areas. We assume an average glacier density of 850 ± 60 kg m -3 (Huss, 2013). According to the glacier inventory RGI 6.0 (RGI Consortium, 2017), the area has a glacierized area of 9,679 km 2 .

Other auxiliary data
To analyse the impact of temperature and precipitation on GS and water mass balance here, we adopt two types of datasets, the gridded reanalysis products and in-situ measurements from four meteorological stations (their locations are 150 labelled in Figure 1, and coordinates are listed in Table S1). Precipitation and temperature records for each site from 2003 to 2016 ( Figure S4) are available from the China Meteorological Data Service Center (http://data.cma.cn/data/weatherBk.html).
Only four in-situ temperature records may not represent the overall condition of the glacierized zone, so we adopt the gridded temperature product from the ERA5 reanalysis data processed by the European Centre for Medium-Range Weather Forecasts (ECMWF). The data is available at https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5. The gridded data is 155 compared with station observations and the correlation index ranges from 0.69 to 0.82 in the interannual variation ( Figure S6), indicating a good consistency. The average values in the glacierized zone from the ERA5 temperature product will be used to represent the temperature condition here.
Global gridded precipitation data Tropical Rainfall Measuring Mission (TRMM) (Huffman et al., 2014) is used to examine the influence of precipitation on water storage. The data is available at https://pmm.nasa.gov/data-160 access/downloads/trmm. Although such a global product is unable to capture the localized spring precipitation in our study area (shown later), it can be used for the investigation of large-scale monsoon precipitation. We also use High Asia Refined analysis (HAR) precipitation product generated using the atmospheric model WRF (Maussion et al., 2014). In this product, long-term precipitation trends are not recommended, but its 10 km spatial resolved seasonal variation is informative to investigate the spatial extent of spring precipitation. The HAR data is available at http://www.klima-ds.tu-berlin.de/har/. 165 Moderate-resolution imaging spectroradiometer (MODIS) data MOD10 (Hall et al., 2006) is used to investigate snow coverage here. The MOD10CM product has a temporal resolution of 1 month and spatial resolution of 0.05-degree. The land surface model Global Land Data Assimilation System (GLDAS)-NOAH (Rodell et al., 2004) is adopted to inspect soil moisture changes, which can be compared to changes in total terrestrial water storage estimated by GRACE. Here, the version 2.1 monthly product with 1.0-degree spatial resolution is used (available at 170 https://hydro1.gesdisc.eosdis.nasa.gov/data/GLDAS/GLDAS_NOAH10_M.2.1/).

Spring precipitation and mass increase
The method of this study is based on the fact that the change in GS mass driven by spring precipitation precedes the change in hydrological signals. Therefore, before introducing the method, we want to demonstrate that GRACE can detect mass changes caused by spring precipitation. At two out of four stations (Bomi and Chayu), spring precipitation is noticeable, 175 even surpassing the summer/autumn precipitation brought by the Indian monsoon ( Figure S4). Yang et al. (2013) provided precipitation records at 22 sites in a broader area and outlined the boundary of the impact zone of the spring precipitation, which roughly covers the glacierized area studied here. Summer precipitation and its associated hydrological mass change are enormous and well recognized, while the spring one is not. Therefore, here we only use the TRMM and HAR results from January to March in Figure 2 to show the initiation of spring precipitation. The precipitation begins to spread south and west 180 since April, when the monsoon gradually increases (not shown here). The TRMM results show a boundary along the latitude 29° N, where the precipitation suddenly decreases to the north. This boundary of change is irrelevant to the terrain and seems to be artificial. This phenomenon cannot be found in the HAR result, which shows abundant precipitation in the glacierized zone in these months. The bottom plots give the GRACE monthly mass anomalies from March to May (two months later than the precipitation), because we found such a time lag in the response of mass change to precipitation. An earlier mass increase 185 from April can be identified in the southeastern part of the Tibetan Plateau.
The performance of TRMM and HAR is compared with our station measurements in Figure S5. According to the in-situ records, the spring precipitation, as a part of the bimodal variation, is obvious at the Bomi and Chayu stations. TRMM is capable of revealing the condition at Chaya at 28.65° N, but performs poorly in regions north of 29° N. The HAR data seems to slightly underestimate the precipitation in April but overestimate in Autumn and Winter, and thus does not present a clear 190 bimodal change.
These results show that spring precipitation can be captured to a limited extent by various reanalysis products and the 'spring-accumulation' pattern of GS mass change in the SETP is recognizable in GRACE observations. The amplitude and phase of the seasonal mass variation from the equivalent water height (EWH) of GRACE are compared in the background of Figure 1. The seasonal amplitude has a spatial distribution similar to that of the Indian monsoon affected area. This pattern 195 reflects the predominance of the monsoon-controlled hydrological process and the weaker glacial signals in this region.
However, the peak month of seasonal changes (the contours in Figure 1) divergently appears earliest in June in the NTM and gradually delays to August in the southern Himalayas, where the annual amplitude reaches its maximum. The shift in peak months reflects the increasing/decreasing contribution from the sinusoid of the hydrological/GS seasonal variation. A key point to point out is that their peaks have a three-month time window offset (we will demonstrate it later), which is a quarter 200 of the annual oscillation cycle and means that the two signals are mathematically orthogonal.

EOF analysis of GRACE
GS and hydrological mass changes dominate the seasonal gravity signals observed by GRACE in this region and they are mathematically orthogonal due to different phases. Therefore, we employ the EOF technique (see the supporting material 205 for mathematic expressions) (Björnsson and Venegas, 1997) to decompose hydrological and glacial signals in the GRACE datasets ( Figure 3). We thus extract two modes with significantly higher explained variances than the other modes (i.e., two significant modes are obtained). Results of different datasets and filters show good consistency, indicating that the first two modes are robust.
Each mode consists of one EOF (the spatial pattern) and one PC (the temporal evolution). Only the first two modes, 210 respectively accounting for 79 ± 5% and 12 ± 4% of the total variance explanation, are shown. Although the first mode is much stronger than the second mode (because the second one is more localized), their signal strength in the glacierized region is comparable on both seasonal and secular temporal scales. Modes above 2 are weak and irregularly show much noise, so they are discarded here.
The trends of the GRACE observation and its decomposed modes are shown in Figure 4. The GRACE observation shows 215 a significant mass loss, which is divided into the first two modes. In the glacierized zone, ~ 2/3 of the negative trend comes from the 2 nd mode and ~1/3 comes from the 1 st mode. The trend of higher modes (> 2) is quite weak (Figure 4d).
According to the spatial coverages (EOF1 and EOF2) and their temporal variations (PC1 and PC2), the first mode covering the low altitude areas on the south of the plateau with a peak month in August/September seemingly represents hydrologic signals and the second mode concentrating in the glacierized region with a peak month in May (the peak month of June in 220 Figure 1 is the mixed result of the first two modes) seemingly represents glaciers. We will verify these hypotheses below.

GS mass estimation from mode 2
GRACE results only show smooth mass patterns and we need some strategy to recover the original amount of mass changes. If we adopt the second mode to estimate GS mass change, this step is necessary. Therefore, a forward modelling method (Yi et al., 2016) is chosen to recover the mass in a pre-defined region iteratively. This method has been widely used 225 (Chen et al., 2015;Wouters et al., 2008), especially in the study of polar ice sheets. In the first step, we divide the glacier mask based on the glacier distribution recorded in RGI 6.0 (RGI Consortium, 2017) (Figure 4e). The lattices have a resolution of 0.5° by 0.5° and are located in glacierized area (by this way we assume the snow signal also comes from the glacierized area, but it does not influence the total mass estimates). In the second step, the mass in each lattice is iteratively adjusted until its smoothing signal (Figure 4f) well matches the GRACE observation (Figure 4c). The details of each combination of datasets 230 and filters are presented in Figure S7 and S8. Therefore, we solve the mass in each combination ( Figure S8). The mass is multiplied by the PC2 series to derive the glacier mass series, and their average is taken as the mass estimate, which will be compared with ICESat observations to test our hypothesis on its physical meaning.

Validation of mode 1 by soil moisture and precipitation datasets
To validate the hypothesis that the first mode represents hydrological signals, we compare it with EOF decomposition 235 results of two other datasets, soil moisture from GLDAS/NOAH and precipitation data from TRMM ( Figure 5). To make them comparable to GRACE in terms of spatial resolution, they are expanded into spherical harmonics, truncated at degree 60, and smoothed by the same filter. Their results are shown in Figure 5. Different from GRACE that has two significant modes, they only have one due to the lack of a glacial signal. The EOF1 of GLDAS/NOAH and TRMM is consistent with that of GRACE.
The PCs are compared at interannual and seasonal scales as well. Note that precipitation is an instantaneous amount, while 240 water storage is a state value, so the former should be integrated in time to make it comparable to the latter. Here, we integrate precipitation in successive six months by an empirical weight function of (1, 2, 3, 4, 5, 6), which will be normalized, and the value is attributed to the sixth month. Different integration methods are tested in the supporting materials.
Of note, mass contributions from the Brahmaputra River and groundwater are absent (and they are troublesome to obtain) and precipitation is assumed as the dominant driver of water storage change without considering the influence of runoff and 245 evaporation (Humphrey et al., 2016), so we do not expect that we can reach a thorough agreement between different datasets. This is acceptable if their temporal consistency is targeted. However, long-term trends in runoff, evaporation and groundwater cannot be ignored and they are differently reflected in these three products, so their trends have been removed before the comparison. The exclusion of unavailable surface water and groundwater in the GLDAS result also causes a weaker strength of its EOF1 compared to that of GRACE. We conclude that these datasets should be comparable in terms of seasonal and 250 interannual variations and the pattern of spatial distribution, but not in the long-term trend and the amplitude of the spatial distribution. The good resemblance in both the EOF1 (spatial pattern) and PC1 (seasonal and interannual temporal evolution) between GRACE, GLDAS/NOAH and TRMM indicates that they reflect similar geophysical processes, i.e., hydrological variations.

Method feasibility and reliability 255
The phase difference of 3 months is a prerequisite for this method and can be verified retrospectively. We tested different phase differences between hydrological and GS signals and decomposed them by the EOF method (refer to section 3.1 in the supporting material). Two conclusions are obtained: only when the GS mass change peaks in May (3 months before the peak month of the hydrological signal) can our simulated result agree with the GRACE observation; the EOF decomposition can well restore both seasonal variation and the trend of the GS signal if the orthogonality is satisfied. 260 Only hydrological and GS signals can explain the first two modes considering their spatial and temporal patterns.
Atmosphere contribution has already been removed in GRACE observations (Dobslaw et al., 2017) and mass transports of solid earth are unlikely to have such strong seasonal variations. We cannot quantify the contribution of groundwater in the second mode, but groundwater is apt to be modulated by stronger rainfalls in summer (Andermann et al., 2012), rather than snowfalls in winter-spring, and groundwater activity will be reduced in winter-spring when the ground is frozen. Therefore, 265 the groundwater component is inclined to exist only in the first mode. The negative trend in the first mode is mostly due to decreasing precipitation in recent years ( Figure S9) and intense groundwater pumping (Shamsudduha et al., 2012). The negative trend in the second mode is supposed to represent GS melting and can be used for estimating GS mass balance.

Glacier and snow mass balance 270
The glacier surface elevation changes measured by the ICESat are compared with our GRACE-based estimates. We interpolate the series of GRACE estimates (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017) into the observation epochs of ICESat (2003ICESat ( -2009 and plot mass changes by GRACE as a function of elevation changes by ICESat (Figure 6a (Figure 6b), while the GRACE PC2 series show a moderate increase after October. We speculate that the snow height does not increase much in the first few months so the contribution of snow mass is not significant. The observations in March and June, as expected, are well above the line, implying an extra snow mass contribution, which can be inferred from the point-to-line vertical distance. The snow contribution relative to the total mass anomalies varies drastically between 0% 280 and 62% with a mean value of 38% within our observation time windows.
The difference between GRACE and ICESat-based estimates of mass change indicates that the snowpack outside the glaciers is a non-negligible contributor to the seasonal mass variation. This is quite different from previous glacier trend estimates, where non-glacier snow was neglected. Based on MODIS observations, the snow coverage area in this region varies from approximately 80,000 km 2 in winter to 30,000 km 2 in summer, both of which are much larger than the inventoried glacier 285 area (Figure 6b). However, heterogeneous snow depths (Das and Sarwade, 2008) and densities across the vast and rugged area make it difficult to measure their mass change by a non-gravimetric way. October with slightly varied initiation and duration from year to year. The maximum mass increase (10-20 Gt) usually occurs in April, when the spring precipitation peaks, and the severest mass loss (-15 --30 Gt) usually occurs in July when the temperature peaks. As the temperature rises from April to July, the monthly mass change curve drops steeply from the peak down to the trough, but the ascending process with mass accumulation is relatively moderate and continuous.
We calculate annual mass increase and decrease by the difference of mass anomalies between November and May and 295 between June and October, respectively. From 2002 to 2017, the annual mass decrease ranged from -43.4 Gt to -73.1 Gt with an average of -57.6 Gt, and the annual mass increase ranged from 35.7 Gt to 63.4 Gt with an average of 50.5 Gt. The seasonal GS mass changes postpone the runoff of ~50 Gt of winter-spring solid precipitation for several months. This amount plays a vital role in the annual streamflow (130.7 Gt on average) of the upper Brahmaputra (Lutz et al., 2014) and is almost ten times the annual net meltwater. Without the buffering effect of the seasonal variation, there will be a tremendous reduction in the 300 streamflow in summer and autumn, when the water demand is high, and adaptive management on the dams in the Brahmaputra will be required to reduce seasonal irregularities in the streamflow (Barnett et al., 2005).

Quantifying the sensitivity of glacier and snow melt to temperature
Temperature is a dominant factor influencing the melting of glaciers (Cogley et al., 2011). Here, the monthly temperature records from the ERA5 product are compared with month-to-month mass changes by GRACE to investigate the sensitivity of 305 the GS mass balance in response to temperature change (Figure 7). Mass changes are negatively correlated with the temperature anomalies by a factor of -1.9 ± 0.2 Gt degree -1 during the ablation season (from May to October) but no correlation is found during the accumulation season (from November to April). The mass peaks around May, when either glacier accumulation or ablation could happen. The temperature averaged in this transitional month is taken as the reference for the temperature anomalies used in the figure and their mass changes are annotated. The highest sensitivity of monthly mass changes in response 310 to temperature is observed in July (3.1 ± 2.5 Gt degree -1 ), when the largest monthly mass loss occurs.
To investigate the impact of climatic variables on the interannual variations of GS mass, we compare annual mass losses (from May to October) with summer temperatures (from June to August) (Figure 7b). The annual mass loss is significantly correlated with the summer temperature, with a slope of -10.7 ± 4.2 Gt degree -1 (P-value: 0.025, R 2 -value: 0.35), indicating that the annual GS mass balance is sensitive to summer temperature. The small value of R 2 is partly due to the relatively large 315 uncertainties of our mass estimate (10 Gt) in this modest range of variation (30 Gt) and the neglect of other factors influencing examined and the SETP shows a widespread high sensitivity with an average value of -1.23 m w.e. degree -1 . Based on the glacierized area of 9,679 km 2 , our estimation is -1.10 ± 0.43 m w.e. degree -1 , which is comparable with the earlier study of Sakai and Fujita (2017). It should be pointed out that annual net mass balance was used in Sakai and Fujita (2017) in 320 comparison with the annual mass loss used in this study, although annual net mass balance is mainly driven by summer melt (Ohmura, 2011).
We could not find a significant relationship between the mass and precipitation changes, probably because our data fail to reflect the strong orographic effect in precipitation, and/or the GS mass gain process is too complex to be attributed to precipitation alone. 325

Comparison with previous estimations on glacier and snow meltwater
The trend of glacier elevation change by ICESat in this study is -0.65 ± 0.20 m w.e. yr -1 during 2003-2009, which lies between the values of -0.30 ± 0.07 m w.e. yr -1 (Gardner et al., 2013) and -1.34 ± 0.29 m w.e. yr -1 (Kääb et al., 2015) in eastern NTM by using alike ICESat dataset (but of an older version), and is close to the trend of -0.62 ± 0.23 m w.e. yr -1 during 2000-2016 by using ASTER (Brun et al., 2017). The trend of GS mass change in this study by using GRACE is -6.5 ± 0.8 Gt yr -1 330 between August 2002 and June 2017. The mass contribution from snow is considerable at the seasonal scale but negligible over 15 years, so the secular trend by using GRACE mainly represents the glacier mass change. Our GRACE trend consists well with the derived glacier mass change of -5.5 ± 2.2 Gt yr -1 by using ASTER (the area-averaged rate in NTM and Bhutan multiplied by the glacierized area of 9,679 km 2 ). In conclusion, both of our ICESat and GRACE estimates agree well with the previous ASTER result in terms of secular trend. 335 A recent result on changes in interannual glacier flow in this region (Dehecq et al., 2018) indicates a strong correlation between ice flow rate and changes in glacial thickness. The interannual variation of GRACE-based mass changes (the 1-year smoothed sequence in Figure 6c) notably shows equilibrium during periods of 2003-2005 and 2011-2014. According to the previous study (Dehecq et al., 2018), thinning glaciers reduce their flow rate by weakening gravitational driving stress; therefore, this balanced mass state may slow down the decreasing flow rate. Coincidentally, we can identify such decelerating 340 phase in the decline of glacier flow rate during 2004-2006 and 2012-2015 (Figure 1 in Dehecq et al., 2018).
GS mass loss is caused by flow, melting, and evaporation processes, while the last one does not contribute to the river flow. Evaporation is important for continental-type glaciers where the climate is usually cold and dry. E.g., it accounts for 12% of the glacier ablation in Tianshan (Ohno et al., 1992). However, the importance of evaporation is greatly reduced in our maritime glaciers due to the extremely humid air and rapid melting. Therefore, we suppose that the mass loss is completely 345 turned into meltwater and can be compared with analogous outputs from models. In our study region, 85% of its meltwater (estimated according to the area proportion) runs into the Brahmaputra and this area accounts for 83% of total glaciers in this basin (9,912 km 2 ). Assuming that the unobserved 17% of glaciers hold the similar rate of GS mass change, our estimate of mass change is scaled by a ratio of 1 × 0.85/0.83 = 1.02 to represent the GS mass change of the entire Brahmaputra Basin.
of Lutz et al. (2014), which showed that GS melt constitutes 33% of the total discharge in the Brahmaputra and that 50% of the annual melt occurs in the summer (Figure 8). GRACE only detects the net change in GS and cannot separate mass ablation and accumulation (see the inset in Figure 8). Because these two processes concur simultaneously in transitional seasons and offset to some extent, the annual mass decrease (total mass loss in a year; here, ranges from 44.5 km 3 to 74.8 km 3 with an average of 59.0 km 3 ) is smaller than the real GS melt. As a result, the annual mass decrease provides a lower bound on annual 355 GS melt each year, rather than an accurate estimate. Instead, the amount of GS melt can be better determined during the summer (from June to August), when the accumulation is supposed to be small. This value can be used to validate the model output. Our result shows that the summer melt ranges from 37.3 km 3 to 62.9 km 3 with an average of 51.6 km 3 , which is over 100% larger than the 23 km 3 GS mass change given in the model of Lutz et al. (2014) (Figure 8). Although extrapolated mass changes for the undetected 17% of glaciers and the neglected summer evaporation may reduce our estimates of summer 360 meltwater, they definitely cannot explain the difference of more than 100%. Among all model estimates, the model of Lutz et al. (2014) reported one of the largest proportions of GS melt contribution (33%), but still largely underestimated the amount of summer meltwater, according to our estimate from satellite observations. Our annual mass decrease (average 49.0 Gt) is still much smaller than the 137 Gt annual meltwater given by Huss et al. (2017). However, this enormous value even exceeds the annual streamflow of 130.7 km 3 in the upper Brahmaputra where all 365 GS meltwater is included (Lutz et al., 2014). The upper streamflow at the Nuxia station (ahead of the main glacier supply area) is ~ 60 km 3 . Therefore, the difference in streamflow between the main glacier supply area is ~ 70 km 3 , and the annual meltwater is unlikely to exceed this value, considering the additional contribution of precipitation. These values generally represent decadal averages at the beginning of this century (Table 1) and they are therefore comparable.

Conclusion 370
In this study, we use GRACE gravimetry to estimate the GS mass balance in the SETP from August 2002 to June 2017.
The second EOF mode of GRACE observations is attributed to changes in GS mass, which can be validated in the following three steps. First, simulation experiment shows that two signals with peaks in August and May can be decomposed unbiasedly by EOF. Second, the first decomposed mode shows consistent spatio-temporal patterns with the soil moisture and precipitation variations from the GLDAS and TRMM data and thus can be reasonably attributed to hydrological processes. Thirdly, the 375 second mode of GRACE signal with a peak in May temporally corresponds to the glacier/snow accumulation and ablation processes and spatially coincides with the glacier distribution, which is also supported by the spring precipitation pattern observed by meteorological stations. Glacier mass change measured by ICESat is further adopted to compare with our GRACE-based GS estimates, and good agreement is reached in the ablation season when the snow contribution is negligible.
The ICESat measurements also show that the seasonal glacier mass variation is large, which is consistent with our finding that 380 GS mass change in this region peaks in May.
The GRACE-based GS mass balance not only shows a long-term decreasing trend of -6.5 ± 0.8 Gt yr -1 , generally comparable with previous studies on glacier mass balance in the SETP, but also newly reveals a strong seasonal variation which postpones water supply of about 50 Gt from winter and spring to summer and autumn. The high sensitivity of glacier mass changes responding to temperature shows that warming climate will exert strong impacts on the glacier and snow mass 385 balance from two aspects. On the one hand, under the current glacier condition, the increase in summer temperature will enhance the annual meltwater by a factor of -10.7 ± 4.2 Gt ℃ -1 ; On the other hand, the seasonal meltwater will shift earlier and reduce its supply in summer and autumn, which is potentially ten times the amount of annual glacier melting. Our estimates of monthly GS meltwater can also give an elaborate calibration on the glacier accumulation and ablation processes in hydrological and glaciological models of the Brahmaputra Basin, which were barely calibrated by GS mass observations and 390 diverged largely on the proportion of seasonal meltwater contribution. Given the high vulnerability to warming temperature, the greater contribution of meltwater to the Brahmaputra streamflow than most model estimates indicates that its water resource allocation will face ominous tension in the future.

Data availability
The data that support this study are mostly publicly open and their sources are indicated in the data section. The 395 meteorological data and the series of glacier mass balance estimate are available upon request to the corresponding author. Figure 1. Geographic environment of the upper Brahmaputra Basin. The boundary of the basin is outlined by the black dashed line. The violet areas in the plateau represent mountain glaciers, but only the darker ones (9,679 km 2 in total) are studied here. The background color shows the amplitudes of annual variation in terms of equivalent water height from GRACE, and their peak months (the month with the peak value in a year) are indicated using contours (e.g. 9 means September). The red triangles mark the location of four meteorological stations. The colored arrows illustrate major climatic 530 factors influencing this region (M: Indian Monsoon; W: Westerly winds; V: Bay of Bengal Vortex). The red box in the inset map marks the location of the study area.