Rapid thinning at a glacier terminus and icefield of 2 Mt . Tanggula , Tibet , detected by TOPEX / Poseidon , 3 and Jason-2 /-3 altimeters : confirmation by satellite 4 imagery 5

An oceanic radar altimeter such as TOPEX/Poseidon (T/P) is typically for observing 15 elevation changes over the open oceans or large inland lakes/rivers, with limited applications over 16 solid earth due to its large footprint and susceptibility to waveform contamination and slope 17 effect. Here we demonstrate that it is possible to construct a long-term time series of glacier 18 elevation change from T/P-series radar altimeters over two flat surfaces near a glacier terminus 19 and an icefield (Sites A and B, with slopes of 2° and 0.8°) in Mt. Tanggula, Tibet, at elevations over 20 5400 m. We retracked radar waveforms using the subwaveform threshold algorithm, selected 21 quality altimeter data (1/4 of the original) with nearly the same slope and adjusted the original 22 elevations by fitting with a time-varying, 2nd order surface. The glacier elevation changes at the 23 two sites from T/P (1993-2002) show seasonal elevation oscillations with linear rates at about -3 24 m/year and abnormal seasonal changes around the 1997-98 El Niño. Site A is over a deep valley in 25 southern Tanggula. Its elevation dropped about 30 m over 1993-2002 (from T/P) and the glacier 26 almost disappeared by 2016 (from altimeters and satellite images). Despite the sporadic Jason-2 27 and Jason-3 altimeter data, we also derived long-term rates of glacier elevation change over 199328 2017. Landsat-derived glacier area and elevation changes near the two sites confirm the rapid 29 glacier thinning from the altimeters. The glacier meltwater near Site A supplied increasing source 30 water to Chibuzhang Co west of Mt. Tanggula, contributing partially to its accelerated rising lake 31 level. The glacier thinning at Site B (icefield) was correlated with the increased discharge of the 32 Tuotuo River in eastern Mt. Tanggula, a source region of the Yangtze River. The successful 33 detection of glacier thinning at the two sites shows that T/P-series altimeters can serve as a virtual 34 station at a flat glacier spot to monitor long-term glacier elevation changes in connection to climate 35 change. This virtual station concept is particularly useful for inaccessible glaciers, but its 36 implementation faces two challenging issues: increasing the volume of quality altimeter data and 37 improving the ranging accuracy over a targeted mountain glacier spot. 38


Introduction
Tanggula Mountains (Mt.Tanggula), at an average elevation of more than 5400m above sea level, are located at the borders of Tibet and Qinghai, largely covered with glaciers, seasonal snow and permafrost.Mt.Tanggula is the host to summer-accumulating mountain glaciers, and its meltwater flows into neighboring lakes such as Chibuzhang Co and Dorsoidong Co [1,2].The glacier meltwater of Mt.Tanggula also contributes to the source water of the Yangtze River, the largest river in China [3].Recent rises of temperature around Mt. Tanggula have accelerated melting of glacier here and raised the concern of lost water storage capacity in Yangtze's source regions [3,4].The melting of Tanggula glaciers is a typical example of worldwide glacier melts caused by global warming.Glacier melts can raise sea level [5], modify ecosystems and alter downstream hydrological systems, among other effects [6].
There have been several studies on the glacier state of Mt.Tanggula using satellite imagery.
For example, Pu et al. [7] examined the mass balance of the Dongkemadi glacier located near Mt.
Tanggula, concluding that the trend of mass balance was positive during the period of 1989-1993, and was negative from 1994 onward, except for the large positive snow accumulation in 1997.In particular, the Dongkemadi glacier experienced a sharp decline since 1998, which was attributed to the increased glacier melt in summer.Furthermore, Qiao [1] used satellite imagery, digital elevation models (DEMs), and glacier inventories to analyze glacier area changes around Dongkemadi over different time spans, showing that the glacier area over the Dongkemadi glacier had declined persistently since the beginning of their data records.
Satellite imagery is useful for detecting glacier area change, but it is limited in providing volume change when no elevation information is available.Satellite altimeter is one of many remote-sensing sensors that support the opposite function of satellite imagery: it can detect glacier elevation change, but not volume change without area information.There are two kinds of satellite altimeters, one based on radar and the other based on laser.There have been numerous radar-based satellite altimeters since the Seasat mission of 1978; information about these missions can be largely found on the popular altimeter web pages of the European Space Agency (ESA; https://www.esa.int/),Archiving, Validation and Interpretation of Satellite Oceanographic data (AVISO; https://www.aviso.altimetry.fr/) and the Jet Propulsion Laboratory (JPL; https://www.jpl.nasa.gov/).In contrast, there are only two laser-based altimeters to date, which are the Ice, Cloud, and Land Elevation Satellite (ICESat) mission (https://www.nasa.gov/mission_pages/icesat)and its successor ICESat-2 (launched in September 2018).Both radar-and laser-based altimeters have been used to measure elevations over ice-covered surfaces.Because of its small footprint (about 70 m), ICESat has been used in several studies of mountain glacier elevation changes [4,[8][9][10].In comparison to a laser altimeter, a radar-based altimeter has much larger footprints, depending on the surface roughness.Because of the large footprints, radar altimeters are largely used over ice sheets, which are sufficiently flat for an altimeter to measure precise elevations.For example, Wingham et al. [11] measured elevation changes over Antarctic glaciers from 1995 to 2006 using the ERS-2 and Envisat radar altimeters.By contrast, only few radar altimeter studies of non-polar glaciers exist.One such example is Lee et al. [12], who used TOPEX/Poseidon (T/P) and Envisat altimeters to detect glacier elevation changes over the Bering Glacier system in Alaska.However, the glacier surfaces studied by Lee et al. [12] were flat and with favorable conditions for precise elevation measurements by the T/P and Envisat altimeters.In addition, over flat croplands in California's Central Valley, North China Plain, and central Taiwan, radar altimeters have been used to determine land subsidence rates at the cm/year accuracy [13,14].
Compared to a radar altimeter like T/P that has an exact repeat period, ICESat's repeat periods are not exact over Tibet.Thus, ICESat and T/P measure glacier elevation changes in different ways.This is illustrated in Figure 1.At an earlier epoch (marked in black), the three elevation differences (∆ℎ , ∆ℎ , ∆ℎ ) are the differences between the ICESat-measured elevations and the elevations derived from a DEM such as one from the Shuttle Radar Topography Mission (SRTM; http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp, see Section 3.1).The mean of the three elevation differences is the representative difference (or glacier elevation anomaly) for this epoch.
At a later epoch (marked in red), the glacier surface has declined, resulting in another three elevation differences (∆ℎ , ∆ℎ , ∆ℎ ) at different locations and a new representative difference.
ICESat uses such representative differences to estimate the rate of glacier elevation change.Figure 1b shows the ground tracks of T/P at an earlier epoch (in black) and a later epoch (in red).The two T/P tracks can have a lateral offset of up to 1 km (from P1 to P2), which is much less than the crosstrack distances of ICESat (Figure 1a).Although P2 is not exactly identical to P1 (the measured elevation at P1 is ℎ ), the measured elevation at P2 (ℎ ) can be easily reduced to an elevation at P1 (P2') using the terrain slope around here.Such slope-corrected elevation measurements (ℎ , ℎ ⋯) can be used to estimate the rate of glacier elevation change.Unlike the T/P tracks (Figure 1b), the ground tracks of ICESat (Figure 1a) are not repeated and the ICESat-derived elevation changes rely on both the ICESat measurements and a DEM.Thus, errors in the DEM, which may exceed 10 m in mountainous area [15], could propagate to errors in the ICESat-derived elevation changes.
In theory, a repeat radar altimeter series such as T/P, Jason-1 (J1), Jason-2 (J2) and Jason-3 (J3) can provide along-track, long-term (1992 to present) glacier elevation changes at mountain glacier sites where elevations can be measured precisely and repeatedly.However, so far there is only one such study [12] over mountain glaciers.The rareness in radar altimetry study of mountain glacier gives rise to two questions: (1) can a radar altimeter like T/P detect elevation changes over alpine glaciers such as those over Mt.Tanggula?(2) if T/P can, how should its data be selected and processed to construct a reliable time series of glacier elevation change?
As it happens, a pass of T/P #155 travels through glaciers in southeastern Mt.Tanggula.
Another pass #242 flies over Chibuzhang Co, whose source water is at Mt. Tanggula (see Section 4.1).The state of Mt.Tanggula glaciers may have affected the state of Chibuzhang Co lake level [3]; the latter can be assessed by the altimeter measurements from pass 242.Thus, the T/P altimeter data from these two passes are ideal for carrying out experiments to answer these two questions and to assess the consequence of Mt.Tanggula glacier change.This paper will also show the links between the rapid glacier decline of Mt.Tanggula, the rising lake level of Chibuzhang Co and the increased river discharge in the source region of the Yangtze River.Because there is no ground truth and no crossover point to validate glacier elevation changes from the T/P-series altimeters near Mt.
Tanggula, we will use changes in glacier elevation and glacier area from Landsat images to assess the altimeter-derived result.

Radar Altimeter Data from the T/P-Series Satellites
In this paper, we used the altimeter data along passes 155 and 242 of the T/P-series altimeters.
The data were downloaded from the AVISO website (see Section 1).The original elevations from AVISO were improved using the waveform retracking method in this paper.A short introduction to the altimeter data we used and the waveform formats is given below.The T/P mission, launched in August 1992, is the first of the T/P-series satellites.We used T/P data only from cycles 11 to 364 in this study.The T/P data consist of geophysical data records (GDRs) at 10 Hz (~660 m along-track sampling interval) and sensor Data Records (SDRs) with 64 waveform gates.The retracking gate is 24.5 for T/P.We also attempted to use data from the J1 mission, which was launched in September 2002 and is the second of the T/P series.However, only few Jason-1 data over Chibuzhang Co can be used in this study, probably because the pre-set heights around Mt. Tanggula were not within the tracking windows of J1's radar altimeter.In addition, several problems in the data processing of J1 rendered unreliable data over land [13][14][15][16].In this study, the J1 altimeter data were used to fill lake level data gaps over Chibuzhang Co between the T/P and J2 records.
The J2 mission was launched in July 2008 and is the third satellite in the series.We used J2 data from cycles 1 to 320 in this study.Like J1, J2 data consist of sensor geophysical data records (SGDRs) at 20 Hz (~330 m along-track sampling interval).A waveform of J2 has 104 gates and its default tracking gate is 32.5.Finally, the J3 mission was launched in January 2016 and is the fourth of the T/P-series satellites.We used J3 data from cycles 1 to 51 in this study.The waveform format of J3 is the same as that of J2.

Satellite Images from The Landsat-5/7 and Sentinel missions
In this study, the satellite images for estimating glacier area changes are from Landsat-5 and -7.
The Landsat-series missions provide satellite images at a 16-day interval.Landsat-

Detecting Glacier Elevation Change by T/P-Series Altimeters
As stated in Section 1, pass 155 of the T/P-series satellites is over Mt.Tanggula.In theory, this pass can provide repeated, along-track glacier elevation measurements every 10 days over many glacier spots around Mt. Tanggula.But the reality is that, we can obtain reliable glacier elevation changes only at glacier spots where the slopes are sufficiently small and the selected elevation measurements can meet the conditions described below.In this paper, we identified two such spots, called site A and site B, to determine glacier elevation changes (Table 1).Figure 2 shows the locations of these two sites, the terrain and three catchment areas that provide source water to lakes and rivers around Mt. Tanggula.The terrain in Figure 2 is based on the Shuttle Radar Topography Mission (SRTM) DEM (in the EGM96 orthometric height system).SRTM DEMs are available at grid resolutions ranging from 1" to 15'.In Tibet, the finest grid resolution is 3" in version 2 prior to 2015.
From our visual inspections of optical satellite images on Google Earth, Landsat images and the contours in Figure 2c and 2d, we conclude that Site A is close to a glacier terminus and site B is over an icefield.Table 1 shows the geodetic coordinates and the orthometric heights of the reference points for the two sites (see below for the definition of a reference point) and the elevation of Chibuzhang Co.In addition, the mean along-track terrain slopes at Sites A and B are 2.3° and 0.8°, respectively.Such small slopes make possible a successful construction of a time series of glacier elevation change at each of the two sites.A glacier height (in the orthometric height system, compatible with that of SRTM DEM) from altimeter measurements can be expressed as where  is the glacier elevation, ℎ is satellite's geodetic height,  is the radar range measurement,  is the geoidal undulation and  is a term containing the following corrections: The first six terms on the right-hand side of Equation ( 2) are corrections for center of gravity, wet tropospheric delay, dry tropospheric delay, ionosphere, solid earth tide and pole tide, which are provided in the GDRs of T/P and the SGDRs of the Jason series of satellites.The last three corrections are the radar range correction by waveform retracking ( ), the slope correction ( ) and height reduction for the geographical distribution of observation points ( ), which are determined in this paper and are explained in the methods below.(In fact, because a major uncertainty with  , this correction is not used, see below) We used two methods to determine glacier elevation changes at sites A and B from the T/Pseries satellites.Method 1, called the standard processing method, uses the approach of Lee et al. [12].Method 2 is called the improved method, in which altimeter-measured elevations over glaciers are improved by waveform retracking, data selection and other extra steps (see below). Figure 3 shows the procedures for both methods for processing the T/P-series altimeter data described in Section 2.1.In Figure 3, the algorithms for "surface gradient correction" and "smoothing" for Methods 1 and 2 are the same.However, Method 2 adds few extra steps to select quality elevation measurements.

Method 1 Method 2
It turns out the slope corrections over Mt.Tanggula are too large to be reliable.Thus, both Methods 1 and 2 do not apply slope correction.This is explained below.First, the slope correction, , is [20 where  is explained in Equation ( 1), and  is the slope of the glacier surface.Around Sites A and B, we have determined the slopes around the radar measurement points using the 3" SRTM DEM (Table 1 and Figure 2).We found that the  values at the raw radar measurement points (for all repeat cycles) can easily exceed 1000 m around Sites A or B. Other rugged terrains around Mt.
Tanggula glaciers can have terrain slopes much larger than 2°, which result in unrealistically large  values.Because  may create huge uncertainties in the radar-measured glacier elevations, we decided not to apply this correction in this paper.Instead, we carefully select elevation measurements that have the same slopes around the reference points of Sites A and B. Such elevation measurements have nearly the same slopes and nearly the same slope corrections (at a site), thus their slope effects can be reduced or even cancelled when differencing such elevation measurements for glacier elevation changes.This reduction of slope effect is an assumption that can be challenged but is useful for determining elevation change from repeat radar altimeter measurements with nearly the same slope over a mountain glacier spot.
The height reduction  assumes that the difference between the elevation at a radar measurement point and that at the reference point (Table 1) is equal to the difference between the corresponding DEM-derived elevations [12][13][14][15][16][17][18][19][20][21][22].Under this assumption, the height reduction is (for both Methods 1 and 2): where ( ,  ) and (, ) are the elevations at the reference point (at latitude  and longitude  ) and the measurement point (at , ) from a DEM, respectively, directly provided in the GDRs and SGDRs without considering the off-nadir effect to the coordinates.In this paper, we chose to use the 3˝×3˝ (90 m) version of the SRTM DEM for the height adjustment in Equation ( 4).
The needed elevations were computed from this SRTM DEM by the bilinear interpolation.Over high mountain glaciers like those in Mt.Tanggula, snowfall may affect the surface roughness and a radar's penetration depth, which will in turn affect the precision of elevation measurements from the T/P-series satellites [12,14,16,21,23].Hwang et al. [14] and Frappart et al. [16] pointed out that Ku-band backscattering coefficients from T/P may be used to correct for the backscattering effect on radar-measured elevations.However, Lee et al. [12] showed that there was no strong correlation between glacier elevation changes and T/P-series satellites' Ku-band backscattering coefficients.As such, this study did not consider the effects of elevation change caused by variations in backscattering coefficient.The last step of data processing in Figure 3 is smoothing.Both Methods 1 and 2 use the Gaussian filter with a window size of 0.5 years to smooth the time series of glacier elevation change from the T/P-series altimeters.
As shown in the result section (Section 4.1), the glacier elevation changes from the two methods (Figure 3) are substantially different.The causes of the differences are attributed to the two different steps marked in blue (retracking) and red (quality data selection) in the methods.The causes are explained below.

(a) Waveform Retracking
Method 1 uses the threshold algorithm [24] to improve the range precision of the T/P-series altimeters.This algorithm is based on the OCOG method [25], which is a purely statistical method.
Lee [22] suggested that a threshold value of 50% should be used over ice surfaces in the OCOG method.In contrast, Method 2 uses the subwaveform threshold algorithm [26] for waveform retracking.The subwaveform algorithm first determines a subwaveform of the full waveform that has the maximum correlation with a reference waveform [27].The selected subwaveform is then used to determine the retracking gate by the threshold algorithm [24].We experimented with three threshold values, 10%, 20% and 50%, in the threshold algorithm.The use of a 50% threshold value is consistent with the suggestion of Lee [22] and results in glacier elevation changes that have minimum fluctuations, clear seasonal signals and trends of glacier elevation change correlated with the trends of glacier area change (Section 4.2).Thus 50% is the adopted threshold value for our final altimeter result.More details about the subwaveform threshold algorithm can be found in Yang et al. [26].

(b) Quality data selection
The quality data selection in Method 1 is based on a simple 3-sigma criterion to remove outliers in the raw radar elevation measurements.That is, the standard deviation of the radarmeasured glacier elevations and their mean elevation for each cycle in a bin around the reference point (Table 1) are computed.An elevation measurement is rejected if its residual (raw elevation minus mean) exceeds three times of the standard deviation.A new standard deviation and mean are then determined without the outliers to detect more outliers.The detection is iterated until no outliers are found.
In comparison to Method 1, Method 2 selects quality data in three steps: (1) removing erroneous elevation measurements, (2) selecting only elevation measurements on the same slope as that at the reference point (Table 1) (3) fitting the selected measurements by a 2 nd order surface while removing outliers.The three steps are explained as follows. Step

1: removing erroneous measurements
This step is based on the recommendations of Kääb et al. [8] and Chao et al. [4], who compared altimeter-measured elevations with those from a SRTM DEM to identify erroneous elevation measurements.Following their recommendations, an elevation measurement is considered erroneous and is removed if its absolute difference with the 3" SRTM DEM-derived elevation exceeds 150 m.
Step 2: selecting only elevation measurement on the same slope In this step, we selected radar elevation measurements that are on the same slope, so that the slope effect can be canceled (see Equation ( 3)).In addition, glaciers not on the same slope may not have the same rate of elevation change [28] and they should not be used for a time series that is assumed to have a uniform changing rate.As shown in Figure 2c and 2d, pass 155 of T/P is not parallel to the directions of the steepest descents of the glacier terrains around Sites A and B. It is known that T/P's repeat tracks can be off by 1 km with respect to a mean track.Thus, a raw elevation measurement can be over a spot where the slope (and slope correction) is substantially different from that near the reference point (Table 1), making the elevation measurement unsuitable for constructing an elevation time series.We experimented with several methods to select sameslope elevations.We reached a simple method for this as follows.Let be the integral number of the elevation of the reference point divided by 100 m.At a radar then the altimeter elevation measurement is selected for further processing.In Equations ( 5) and ( 6),

𝐻
is the elevation of the reference point at Site A or B (Table 1) and  is the elevation interpolated from the 3˝×3˝ SRTM DEM at point i.In this step, each of the selected raw elevations (after step 2) from all repeat cycles in a bin is fitted by a function in time and space in the following observation equation [29,30]: where  (, , ) is glacier elevation from repeat cycle j at point i with geodetic latitude ( ), longitude () (same as those in Equation ( 4)), and measurement time (),  is the residual,  and  are the latitude and longitude of the reference point in a bin (Table 1),  is the mean elevation of the bin,  is the initial rate of elevation change, and    are the coefficients of a 2 nd order surface model.In Equation ( 7), we assume that the altimeter elevation measurements from all repeat cycles must fall into a 2 nd order surface that linearly changes with time.The 7 coefficients in Equation ( 7) were estimated by the method of least-squares requiring that the weighted sum of the squared residuals be a minimum.The weight for a raw measurement  (, , ) is the inverse distance of this point to the reference point (Table 1).The actual least-squares method we used is robust and iterative.That is, after initially least-squares estimating the 7 coefficients in Equation ( 7), for every measurement we examined whether  is three times larger than the a posteriori standard deviation.If this happened, the measurement was considered anomalous and was removed.After all such anomalous measurements were removed, we estimated the coefficients again.This process of parameter estimation and outlier removal was iterated until no outliers were found.After the surface fitting, the adjusted elevation for  (, , ) was computed as All the adjusted elevations in a repeat cycle j were then reduced to the elevations at the reference point using the height reduction in Equation ( 4): ( ,  , ) =  (, , ) +  (9) By averaging all such elevations, we computed a representative elevation for the repeat cycle j: where n is the number of measurements.The elevation ( ,  , ) was used to construct the time series of glacier elevation change in this paper.

Detecting Lake Level Change by T/P-Series Altimeters
Pass 242 of T/P travels through Chibuzhang Co, where we computed lake level changes to see the hydrological consequence of Mt.Tanggula's glacier melt (from pass 155) in Basins I and III (Figure 2).Our data processing procedure for lake level change is the same as the one used in Hwang et al. [2], and is summarized in four steps below.
Step 1: Selecting altimeter data over Chibuzhang Co First, the GMT-defined lake polygon was used to determine Chibuzhang Co's coverage and then a window was set for data selection.According to the experience in Hwang et al. [2], the altimeter data for a best result of lake level change should be those located as close as possible to the center of the lake to avoid land interference on altimeter waveforms.Thus, we selected only altimeter data that are 2 km away from Chibuzhang Co's shores.
Step 2: Determining lake elevations from T/P-series altimeters Like glaciers, a lake's ellipsoidal height at a spot is the difference between the ellipsoidal height of the altimeter and the altimeter's range measurement.To improve the ranging accuracy over Chibuzhang Co, we also used the same subwaveform threshold algorithm (see Method 2, Figure 3; Section 3.1) to determine the range corrections for the selected elevation measurements.However, a 20% threshold value was used over Chibuzhang Co, instead of 50% that was used over glaciers.The use of 20% is based on the result of Hwang et al. [2] over many Tibetan lakes.
Step 3: Excluding outliers using the 3-sigma criterion To obtain reliable lake elevation measurements, we applied the 3-sigma criterion to exclude outliers.First, the standard deviation and the simple mean of the selected elevations in Step 2 were computed.We then examined the residual of each selected elevation.If the residual exceeded three times of the standard deviation, the elevation was disregarded and a new mean and a new standard deviation were computed.This iteration stopped when outlier was no longer found.
Step 4: Averaging lake elevations For each repeat cycle, the selected, outlier-free lake elevation measurements were used to compute an elevation at the center of the lake (its position is the same for all cycles).Because the length of pass 242 over Chibuzhang Co is short (Figure 2b) and the selected altimeter data are near the lake center, we neglected geoidal differences when averaging the elevation measurements to obtain the mean elevation.

Detecting Glacier Area Change and Elevation Change by Landsat Imagery
For comparison with the altimeter-derived glacier elevation changes, we determined area and elevation changes around Sites A and B using superimposed Landsat images.First, the two Sentinel-1A SAR images were used to construct a DEM using the procedure given in Figure 4.
From each of the two Interferometric Wide Mode images, we selected one sub-swath image and then aligned the resulting images.Next, we performed the interferometry process to obtain phase differences, debursted the images and applied the Goldstein filtering to the interferograms.Then, we used the SNAPHU software [31][32][33] to unwrap the phases.The unwrapped phases were converted to elevations, which were then corrected for the range-Doppler terrain effect to obtain the final DEM.We determine whether a Landsat image pixel is over glacier using the steps below.
Step 1：Computing the Top-of-Atmosphere (ToA) reflectance We selected cloud-free images over Mt.Tanggula to compute pixel-wise glacier coverages.For each pixel, first we computed the Top-of-Atmosphere (ToA) reflectance using [34]  =  ×  +  (11) where  is the radiance,  and  are the gain and bias for the Green or the mid-infrared band, and  is the digital number of the pixel.Then we determined  =  ×  ×   ×  (12) where  is the ToA reflectance (unitless),  is the distance between the Earth and Sun,  is the mean solar exoatmospheric irradiance, and  is the solar zenith angle.

Step 2：Computing Normalized Difference Snow Index
The ToA reflectance was used to compute the Normalized Difference Snow Index (NDSI) where  and  are the reflectances in the green and mid-infrared bands from Equation (12).
Following the recommendation of Cea et al. [35], we decided that if NDSI ≧0.4, then the pixel is over glacier.The resulting image was then compared with the original image to ensure the glacier coverage is correct.In total, 58 Landsat images were used.The determinations of area change and elevation change were carried out in the following two separate steps: (1) Area change As shown in Figure 5c, for each image the glacier edge is detected using NDSI and the edge defines the glacier coverage at the image-acquiring time.We selected the smallest glacier-coverage area during the dry season of a year to construct time series of glacier area change (Section 4.2). (

2) Elevation change
As shown in Figure 5d, we co-register the glacier edge with the DEM from Sentinel-1A.For each image, the glacier elevations at the pixels near the glacier edge are determined and averaged to produce a mean glacier elevation.The mean glacier elevations from the 58 Landsat images form a time series of imagery-derived glacier elevation changes, for comparison with altimeter-derived elevation changes in Section 3.1 images in 2016, overlapped with the glacier edge to approximate the height around the terminus.

Glacier Elevation Changes from T/P-Series Altimeters
First, we compare the glacier elevation changes derived from T/P-series altimeters by Method 1 (Figure 6) and Method 2 (Figure 7). Figure 6 shows that the year-to-year variations in glacier elevation by Method 1 at Sites A and B are up to 500m, which are too large to be reasonable.By contrast, the elevation variations by Method 2 (Figure 7) are much smaller and the glacier elevations at the two sites declined at a steady rate.Therefore, we believe that only Method 2 can produce a reliable rate of glacier elevation change at the two sites.Method 1 uses a more relaxing selection criterion (see Figure 3) that results in more elevation observations, as compared to Method 2. In addition, the number of raw J2 data points at Site A was small and there were only few J2 data points at Site B (not shown in Figure 7b).Because of the problem with Method 1, we process only T/P data.Thus, Figure 6 shows only the T/P result.In the following, we will only discuss the result from Method 2.
Table 2 shows the rates of glacier elevation change by Method 2 and the rates of lake level change (the discussion about lake level change will be presented in Section 5.1).The rates of glacier elevation change from T/P at Sites A and B are close and both show that on average the glacier elevations dropped about 3 m per year during 1993-2002.Table 2 shows that only 22-27% of the raw T/P altimeter data were selected at the two sites.The percentages of the usable J3 data are slightly lower than those of T/P, but J3's data period (early 2016 to late 2017) is too short to draw a conclusion on J3's usable data percentage.Again, only few raw J2 elevation measurements were found near the two sites.The low percentages (about 1/4 of the original) of the usable T/P and J3 are mainly due to the quality data selection procedure given in Section 3.1 (see Method2, Figure 3).In comparison to the low percentages of the usable T/P data over the two glacier sites, the percentage over Chibuzhang Co is much higher (69%, Table 2).
Despite the data volumes with J2 and J3 data at Site A, we estimate the rate from T/P, J2 and J3 over 1993-2017 and the resulting rate is -3.57±0.15m/year, which is consistent with the T/P-only rate (-3.71±0.30;Table 2) over 1993-2002.At Site B, the rate from T/P and J3 is -1.45±0.177b), which require more investigations about the potential causes.One likely cause of the large oscillations in glacier elevation is the rising temperature and increased rain that may have roughened the glacier surfaces; such surfaces may have led to contaminated waveforms beyond repair by the subwaveform retracking (Section 3.1).Altimeter-detected large oscillations over Bering Glacier in Alaska were also reported by Lee et al. [36].

Comparison with the Landsat-Derived Area and elevation Changes
The Method 2 altimeter results at Site A and B are compared with the Landsat-derived area and elevation changes (see Section 3.3).Figure 9 shows two Landsat images in 1986 and 2015 over Mt.Tanggula.These two images show that Site A is near a retreating glacier terminus to its south and Site B is over an icefield (see also Figure 2d) and near a retreating glacier to its northeast.
Several other glacier terminuses on these two images also show signs of glacier area losses.Table 3 shows the glacier area changes over different periods from 1986 to 2015 from the Landsat imagery.
The glacier area changes near Site A in all periods of 1986-2015 were negative, but the area changes near Site B fluctuated between -0.4 km 2 and 0.3 km 2 .In the processing of Landsat-series imagery, we calculated the minimum of the annual snow/ice area from as many cloud-free images as possible.
We assume that this minimum area was not affected by snow cover.However, some satellite images failed to deliver results over Site A due to cloud covers that made it difficult to identify snow coverages using Equation (13).
Because the footprints of the T/P altimeters at Sites A and B are not right at the glacier edges (note that Landsat-based elevation changes are based on glacier edges), we computed elevation changes at two sub-regions of A (A1 and A2) and B (B1 and B2), as shown in Figure 10(a).Figure 10 b, c, and e show the yearly changes in glacier elevation derived from Landsat.Both Landsat-based elevation changes at A and B show declining trends since 1986, but with magnitudes of thinning smaller than those from altimeters (Figure 7).For example, the Landsat-derived thinning rate is -1.20 m/year (A2), compares to -3.71 m/year (from T/P).The Landsat-derived rates at B are smaller than those at A1 and A2 and have larger fluctuations.Furthermore, both the Landsat and altimeterderived rates show larger elevation oscillations at Site A than B. The potential causes of the differences between the rates from satellite altimetry (Figure 7) and Landsat imagery (Figure 10) are: (1) An altimeter-derived rate of glacier elevation change is the rate over footprints with a 1-km radius [37], whereas the Landsat-derived rate is near glacier edges with elevations from glacierdepleted ground to the border of the area under consideration (see Figure 11a.In summary, both the satellite altimetry and imagery result in this paper show thinning of glaciers at Sites A and B, despite different rates of thinning caused by the natures of these two sensors.The Landsat imagery delivers only yearly minimum glacier areas without seasonal area changes (see Section 3.3).In contrast, the T/P-series altimeters can produce clear seasonal elevation changes because of its 10-day repeated observations (subject to the usable data due to the quality data selection).From Figure 7 and 11, we can draw the following conclusions: (1) the glacier elevation changes at Site A from the T/P-series altimeters (Figure 7a) were positively correlated with the Landsat-derived glacier area changes (Figure 11a); (2) At Site B, the glacier elevation changes from T/P (1993-2002) were positively correlated (Figure 7b) with the glacier area changes(Figure 11b), but it is not clear how the correlation behaved after 2003 because there was no J1 and J2 altimeter data and the J3 record was short.In summary, the Landsat result has largely confirmed the glacier elevation declines at Sites A and B detected by the T/P-series altimeters.

Glacier declines and rises of Chibuzhang Co Lake Level and Tuotuo River Discharge
Chibuzhang Co and Dorsoidong Co are located west of Mt.Tanggula (Figure 2; the two lakes are jointly named Chibuzhang Co in this paper).It is known that the levels of two lakes were affected by the glacier melt of Mt.Tanggula [3,38,39].This study also used the T/P-series altimeter data to detect lake level changes of Chibuzhang Co (Section 3.2); the result is shown in Figure 12.
Both the time series of glacier elevation change (Figure 7, T/P only) and lake level change (Figure 12, T/P, J2 and J3) show clear seasonal variations, but the two time series have opposite trends.In general, the glacier elevation of Mt.Tanggula is the highest in late winter (northern hemisphere), and the lowest in late summer.In contrast, the lake level of Chibuzhang Co is the highest in summer and the lowest in winter.The summer lake level highs of Chibuzhang Co are partially caused by Mt.Tanggula's glacier melt and partially caused by rain in Basin I (Figure 2a).
As shown in Figure 2b, the distance between T/P pass 242 (lake center) and a nearest equilibrium line of Mt Tanggula glaciers in Basin I is about 92 km, and the distance between this nearest equilibrium line and Site A is about 16 km.Although Site A belongs to Basin III (Figure 2b), rather than Basin I, Site A's short distance to Basin I implies that the glacier elevation variations observed at Site A (Basin III) should be highly correlated with those in Basin I.That is, the T/Pobserved glacier melts at Site A may signify glacier melts in Basin I that in turn led to the increase of Chibuzhang Co's lake level seen in Figure 12. al. [3] and Jiang et al. [38]).Furthermore, the large lake level fluctuations from T/P (1993-2002) can be explained by the merger of the two lakes.That is, before the two lakes were merged, only Chibuzhang Co received meltwater from Mt. Tanggula, resulting in large lake level fluctuations due to a relatively small lake volume.As the total lake volume increased by the merging, the same amount of meltwater caused less lake level oscillations in the Chibuzhang Co-Dorsoidong Co lake system.To quantify the effect of lake volume change on the lake level fluctuations, we determined the areas of the three lakes, which are 438, 120, and 393 km 2 for Dorsoidong Co, the middle lake (with pass 242, Figure 2b) and Chibuzhang Co, respectively.This implies that the middle lake is about 23% of Chibuzhang Co.The channel from this lake to the main body of Chibuzhang Co may be blocked occasionally to become an independent lake.It is suspected that, during 1992-2002, the channel was blocked and the middle lake was isolated from main Chibuzhang Co.Because of the relatively small area of the middle lake and the potential isolation from Chibuzhang Co, the middle lake has experienced large lake level fluctuations in the T/P era seen in Figure 12.
It is likely that increased rain and temperature-induced glacier melt have contributed to the rises of lake levels of Chibuzhang Co (Figure 12).This is discussed in several studies.For example, Pu et al. [5] and Chao et al. [4] showed that rain increased at rates of 2.714 and 2.285 mm/year over 1988-1999 and 1966-2013, respectively, at the Tuotuo River meteorological station near Mt.
Tanggula.In addition, Chao et al [4] showed an increasing rate of 0.032 °C/year in temperature over 1979-2007 at the Tuotuo River station.
As discussed in Section 1, the mass balance of Xiao Dongkemadi glacier (near Mt.Tanggula) has been studied by Pu et al [5], who found that the elevation of Xiao Dongkemadi's snow equilibrium line increased in 1997, but decreased suddenly in 1998.Pu et al. [7] speculated that the dramatic elevation oscillations of the snow equilibrium line were caused by the 1997-98 El Niño.
The lake level time series in Figure 12 has an anomalously high value in the summer of 1998.This high may be associated with the sudden elevation rise of the snow equilibrium line in the summer of 1998, which implies more meltwater to enter Chibuzhang Co.In addition, Figure 7 shows several anomalous glacier elevation changes around 1997-1998.First, at Site A (Figure 7a) the glacier elevation was especially large in the winter of 1996-1997, and there was no glacier high in the winter of 1997-1998.Instead, the peak glacier elevation at Site A occurred in the spring of 1998 (in a regular spring the glacier should be low).Furthermore, at Site B (Figure 7b) there was a large glacier elevation drop in the spring of 1998, concurring with a sudden increase of the Chibuzhang Co lake level (see Figure 12).The difference between the patterns of glacier elevation change at Sites A and B around 1997-1998 could be due to the fact that Site B is on the windward side of the east Asia monsoon winds in spring and Site A is on the leeward side.The windward side (Site B) may have experienced large fluctuations in temperature and precipitation around 1997-1998 that led to large glacier elevation oscillations, which in turn roughened the glacier surfaces to increase T/P's ranging noise level.These episodic elevation changes in glacier and lake over Mt.Tanggula are more likely to be detected by a series of repeat altimeter missions than by a non-repeat mission like ICESat.
Our final discussion is on river discharge affected by the glacier melt of Mt.Tanggula.In Figure 13, we overlap the glacier elevation changes (Figure 5b) with the annual discharges [4]    In comparison to ICESat, the T/P-series altimeters delivered 96 effective elevations at site A, which show clear seasonal glacier elevation changes (Figure 7a).Although ICESat was able to detect lake level changes over Chibuzhang Co and Dorsoidong Co [3,39], ICESat delivered only 17 effective lake level observations over 2003-2009 and was unable to detect seasonal lake level oscillations in the two lakes.In contrast, the T/P-series altimeters delivered 246 effective observations over Chibuzhang Co from 1993 to 2017.Thus, if a flat glacier spot or a lake is visited by the T/P-series altimeters, elevation changes here at time scales from months to years can be detected.Finally, by our quality data selection criteria we selected ¼ of the raw T/P measurements for our glacier study.Thus our usable T/P data rate is one observation every 40 days, similar to the data rate of Envisat (one observation every 35 days if all Envisat data are usable over mountain glaciers).

Conclusions
This study is the first of this kind to detect glacier elevation changes and long-term trends over a difficult terrain like Mt. Tanggula by oceanic radar altimeters.The implication is that, with an improved data processing technique, we can use radar altimeters to observe seasonal, inter-annual and secular changes over inaccessible glaciers in regions such as the Himalayas, the Pamirs, the Karakoram, and many other regions covered with mountain glaciers.Such glacier elevation changes will be useful in examining the consequences of regional climate changes and help to analyze the trends of glacier source water for regions that depend on glaciers for water supplies.At Sites A and B, the T/P-series altimeters are like virtual stations that provide continuous, long-term glacier elevation changes in connection to climate change, and help to explain the variations in lake level and river discharge.
Obtaining the result reported in the paper is not without effort: every effort should be made to select the best altimeter data, retrack waveforms for better ranging precisions and correct for the terrain effect.About ¼ of the raw T/P elevation measurements at Sites A and B can be used for glacier elevation determination, but only a small amount of J2 data can be used for this purpose.
The usable data volume of J3 is improved over J2, but is still less than that of T/P.T/P's success in yielding sufficiently large usable data volumes at Sites A and B suggests that the responsible agency of a future altimeter mission should consider using an altimeter similar to TOPEX and fine-tuning the altimeter's tracking and retracking algorithms to maximize the number of effective elevations over a difficult terrain like the glaciers of Mt.Tanggula.In summary, increasing the volume of quality altimeter data and improving the radar ranging accuracy remain two challenging issues when implementing the concept of altimeter-based, virtual glacier station.

Figure 1 .
Figure 1.(a) Three non-repeat ICESat tracks and ICESat-derived elevation changes at an earlier epoch (in black) and at a later epoch (in red) over a glacier, (b) two ground tracks of TOPEX/Poseidon (T/P, from two cycles) belonging to the same repeat pass, with the measured elevations at P1 and P2.See the text for the explanations of the symbols in this figure.The elevations and distances are not to scale.

Figure 2 .
Figure 2. (a) The terrain over Mt.Tanggula, Chibuzhang Co and its sister lake Dorsoidong Co.The glaciers, shaded in white, are defined by the Global Land Ice Measurements from Space (GLIMS) database (https://www.glims.org).The lakes (blue-shaded) are defined by GMT and Google Earth.(b) The ground tracks of T/P pass 155 (over Sites A and B) and 242 (over Chibuzhang Co).The red points show Sites A and B in Table 1.(c) The elevation contours and T/P pass 155 at Site A, and (d) Site B. The red polygons represent catchment basins around Mt. Tanggula.Partially, Basin I supplies source water to Chibuzhang Co, Basin II to the Yangtze River, and Basin III to Selin Co in Tibet (at about 31.81°N,89.01°E).

Figure 4 .
Figure 4.The numerical procedure for generating a DEM from two Sentinel-1A SAR images for Landsat derivation of glacier elevation near Sites A and B.

Figure 5 .
Figure 5.The workflow that identifies glacier area and elevation in each Landsat image covering Site A (glacier terminus).(a) The Landsat-5 image without cloud cover in 1986 (for example).(b) NDSI map of a typical glacier (see Eq. 13).Warmer color indicates higher possibility of ice surfaces.(c) The glacier edge (white line) delineated from NDSI.(d) DEM derived from Sentinel-1A SAR m/year over 1993-2017 (no J2 data were used), compared to -3.08±0.20 m/year from T/P only.There was no winter glacier high in 1997-1998 at Site A and there were large glacier elevation fluctuations at Site B in the spring of 1998 (by T/P) and in the summer and fall of 2017 (by J3).A possible link of such abnormal glacier elevation changes to the 1997-1998 El Niño and the lake level of Chibuzhang Co will be discussed in Section 5.1.

Figure 6 .Figure 7 .
Figure 6.Glacier elevation changes from T/P over 1993-2002 by Method 1 at Site A (blue line), same, but at Site B (red line).

Figure 8 .
Figure 8.A 2016 Google Earth image with an added coordinate frame, showing that Site A (red dot) is over a deep valley, where glaciers below the T/P ground track (in green) almost disappeared by 2016.

Figure 9 .
Figure 9. Two sample images of Landsat over Mt.Tanggula.(a) Landsat-5 image in 1986 and (b) Landsat-8 image in 2015.These two images show reductions of glacier areas near sites A and B and other glacier terminuses between 1986 and 2015.

Figure 11a and 11b
Figure 11a and 11b show the time series of annual minimum glacier area around Sites A and B. At Site A, the glacier areas declined almost monotonically, except the slight increases in 1987, 1998-1999 and 2010-2016.At Site B, the pattern of glacier area changes is different from that at Site A: the mean rate of glacier area change at Site B was negative before the spring of 2005 and the mean rate became positive over 2005-2015.However, the mean rate of area change at Site B over 1986-2015 was negative.In addition, the glacier area changes at Site B have larger oscillations than those at Site A.

Figure 10 .Figure 11 .
Figure 10.(a) The locations of sub-areas A1, A2 (for Site A), B1 and B2 (for Site B) where rates of glacier elevation change are determined by Landsat images, (b) (c) Rates at A1 and A2 over Site A, (d) (e) Rates at B1 and B2 over Site B.

Figure 12 .
Figure 12.Time series of lake level change of Chibuzhang Co from T/P-series altimeters over 1993-2017, overlapped with glacier elevation changes at Site A.

Figure13.
Figure13.Annual river discharges from 1990 to 2012 at the Tuotuo river station near the source region of the Yangtze River (Basin II, Figure2a), overlapped with the glacier elevation changes at Site B (Figure7b).

5. 2
Effective Elevation Measurements From T/P and ICESat Over Mt.TanggulaThe T/P-series altimeters are radar-based altimeters originally designed for oceanic studies.In this study, we extend their oceanic function to detecting elevation changes near a glacier terminus and an icefield over Mt.Tanggula.Although these radar altimeters have shortcomings compared to ICESat, they have detected glacier elevation changes over terrains like those at Sites A and B; such results could be achieved by ICESat only if ICESat had a repeat pass there.Here we briefly compare the results of lake level and glacier height monitoring by the T/P-series radar altimeters and by the ICESat laser altimeter.First, ICESat has the advantages of wide spatial coverage, small illuminated footprint, and high ranging accuracy.However, the GLAS laser onboard ICESat had problems in 2003, so that ICESat had no repeat measurements over Tibet during its entire mission.Thus, one of ICESat's weaknesses is its coarse temporal resolution in Tibet.For example, Chao et al.[4] used ICESat to detect glacier elevation of the entire Mt.Tanggula, resulting in only 17 time-lapsed elevations that were used to determine rates of glacier elevation change from 2003 to 2009.In addition, some regions in Tibet may be visited by ICESat in only one season in a year.

Table 1 .
The geodetic coordinates of the reference points for Sites A and B and Chibuzhang Co

Table 2 .
Altimeter-derived rates of glacier and lake elevation changes a Amplitude of annual variation b Rate from T/P for Chibuzhang Co only c Rate from J1, J2 and J3 over 2003-2017 for the Chibuzhang Co-Dorsoidong Co lake system (see Section 5.1)

Table 3 .
Estimated glacier area changes (km 2 ) from Landsat imagery over different periods

Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 15 November 2018 doi:10.20944/preprints201811.0347.v1
of the Tuotuo River collected at the Tuotuo River gauge station (at 34.22°N, 92.44°E) up to 2012; beyond 2012, no river discharge is available to this study.A detailed discussion on the Tuotuo River discharges and ICESat-derived elevation changes of Geladandong glaciers over 2003-2009 has been given by Chao et al. [4].The river discharges in Figure 13 show an increasing trend that is likely due to the increased glacier meltwater from Mt. Tanggula, as explained by Chao et al. [4].However, ICESat provided glacier elevations only between 2003 and 2009.On the other hand, T/P has provided continuous glacier elevation changes from 1993 to 2002 (Figure 13) that were negatively correlated with the river discharges.That is, T/P has filled the missing glacier elevation records before 2002 and helped ICESat to explain the increased river discharges of the Tuotuo River from 1993 to 2009, potentially due to global warming.If we had glacier elevation changes from J1 and J2 (2003-2016), a long-term relation between glacier melt and river discharge could have been established.