Glacier and lake level change from TOPEX-series and Cryosat-2 altimeters in Tanggula: Comparison with satellite imagery

We detect long-term elevation changes near a glacier terminus and an icefield in Tanggula Mountains from altimeter data from the TOPEX/Poseidon (T/P), Jason-1 (J1), Jason-2 (J2), and Jason-3 (J3) altimeters. An altimeter processing technique is developed for the detection, including waveform retracking, quality data selection, and elevation adjustment. The altimeter-observed glacier thinning is confirmed by the direct elevation differences between the digital elevation models from the satellite missions TanDEM-X and SRTM, and by the glacier area losses from Landsat images. At the two glacier sites, the altimeter-derived elevations show seasonal oscillations, with mean rates of -3.71 ± 0.3 and -3.08 ± 0.20 m year -1 over 1993–2002 (T/P only), and -1.88 ± 0.06 and -1.46 ± 0.07 m year -1 over 1993–2020 (T/P-Jason altime-ters), with abnormal changes around the 1997–1998 El Niño. The declining rates sug-gest that the persistent losses of glaciers in Tanggula Mountains can no longer sustain the large thinning rates. The glacier feeds the nearby Chibuzhang Co lake, and lake level and glacier level changes from Cryosat-2 are consistent with those from the T/P-Jason altimeters. The glacier meltwater contributed to the accelerated rise of the glacier-fed Chibuzhang Co, at a rate of 0.15 ± 0.05 m year -1 over 1993–2002, increasing to 0.42 ± 0.01 m year -1 over 2003–2020. The lake levels of Chibuzhang Co experienced notable drops after the 1997–1998 and 2014–2016 El Niño. The Cryosat-2 result shows the altitude effect of glacier change: the higher the glacier, the less it melts. A repeat altimeter can provide time-lapsed elevation measurements as a virtual glacier station to monitor glacier melt caused by climate change.


INTRODUCTION
Satellite altimeters have been widely used in monitoring glacier elevation changes (Flament and Remy 2012) and continental surface waters (Calmant et al. 2008). There are two kinds of satellite altimeters, one based on radar and the other based on laser. There have been numerous radarbased 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/). Nadir-looking radar altimeters like Envisat have detected water level changes over a 200-m wide and steepsided river (Biancamaria et al. 2017), demonstrating the concept of "virtual" station using elevation observations at a fixed spot along a ground track of a repeat altimeter mission. Despite Cryosat-2's long repeat period of 369 days, the satellite's sub-cycles may also provide data for virtual Terr. Atmos. Ocean. Sci., Vol. 32, No. 1, 1-20, February 2021 stations every 30 days. In fact, there have been a number of studies that used Cryosat-2 altimeter data for glacier elevation detection, e.g., Helm et al. (2014), Siegfried et al. (2014), Gourmelen et al. (2017), Ciracì et al. (2018), and Jakob et al. (2020). Extracting useful heights from an exact repeat mission like Topex/Poseidon (T/P) requires different computational algorithms to a mission like Cryosat-2 than those used over lakes and oceans. In comparison to a laser altimeter such as ICESat (https://www.nasa.gov/mis-sion_pages/ICESat) and its follow-on ICESat-2 (launched in September 2018), where the diameter of the laser footprints are between 15 and 70 m, a radar-based altimeter has diameter of several km and lower accuracies. Consequently, radar altimeters are most frequently used over sufficiently flat ice sheets. For example, Wingham et al. (2009) measured elevation changes over Antarctic glaciers from 1995 to 2006 using the ERS-2 and Envisat radar altimeters. Only few radar altimeter studies over non-polar regions exist. One such example is Lee et al. (2013) who used TOPEX/ Poseidon (T/P) and Envisat altimeters to detect glacier elevation changes over the flat Bering Glacier system in Alaska. Moreover, radar altimeters have also been used to determine land subsidence rates at the cm/year accuracy over flat croplands in California's Central Valley, North China Plain, and central Taiwan (Hwang et al. 2016a).
In theory, repeated radar altimeter series such as T/P, Jason-1 (J1), Jason-2 (J2), and Jason-3 (J3) can provide along-track, long-term (1993 to present) elevation sat mountain glaciers (these four satellites are called T/P-Jason altimeters below). However, at present there is only one such study (Lee et al. 2013) 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 precise elevation changes over alpine glaciers such as those over Tanggula Mountains? (2) If T/P can, how should its data be selected and processed to construct a reliable time series of glacier elevation changes?
To answer these questions, we identify two passes of the T/P-Jason satellites over Tanggula Mountains for our experiments. Pass 155 of T/P travels through glaciers in southeastern Tanggula Mountains. Another pass (pass 242) flies over Chibuzhang Co (lake), whose source water is from Tanggula Mountains. It has been suggested that Tanggula Mountains glaciers may have affected the lake level of Chibuzhang Co (Zhang et al. 2011;Tseng et al. 2016;Jiang et al. 2017). However, there is no direct comparison between glacier elevation changes and lake level changes. Lake level changes of Chibuzhang Co have been studied using data from T/P pass 242 (Hwang et al. 2016b;Tseng et al. 2016).
The objective of this paper is to show a technique for computing glacier elevation changes at two smooth glaciers around Tanggula Mountains from the T/P-Jason altimeters, demonstrating the concept of an altimeter-based virtual gla-cier station. Such elevation changes will be compared with those from the Cryosat-2 altimeter and with those between the elevations from the satellite SAR missions SRTM and TanDEM-X (Rizzoli et al. 2017). We will also discuss the relation between glacier elevation changes and lake level changes around Chibuzhang Co.

GLACIER RETREAT OVER TANGGULA MOUNTAINS
Tanggula Mountains, at an average elevation of more than 5400 m above sea level, are located at the borders of Tibet and Qinghai, largely covered with glaciers, seasonal snow and permafrost. Figure 1 shows two Landsat images in 1986 and 2019 over Tanggula Mountains, which indicate clear glacier retreats over several glacier terminuses. Such retreats will lead to glacier area losses and elevation declines. The meltwater of Tanggula Mountains flows into neighbouring lakes such as Chibuzhang Co and Dorsoidong Co (Qiao 2010;Hwang et al. 2016b), and contributes to the source water of the Yangtze River, the largest river in China (Liu et al. 2020). Recent rises in temperature around Tanggula Mountains have accelerated glacier melting in Tibet and raised the concern of lost water storage capacity in Yangtze's source regions (Liu et al. 2020). The melting of Tanggula glaciers is an example of worldwide glacier retreats caused by global warming. Glacier melting can raise sea level (Zemp et al. 2015), modify ecosystems and alter downstream hydrological systems, among other effects (Bliss et al. 2014). In addition, the level of a glacier-fed lake like Chibuzhang Co in Tibet was shown to expand due to increased supply of glacial melt water (Song and Sheng 2016).
In this paper, we identified two such spots, called sites A and B, to determine glacier elevation changes from the T/P-Jason altimeters with a theoretical sampling of 10 days. Table 1 shows the geodetic coordinates and the elevations of the two glacier sites and Chibuzhang Co. Figure 2a shows the topography around the two sites (based on the SRTM DEM; see section 3.3), and three catchment basins (I to III) that provide source water to lakes and rivers around Tanggula Mountains. The level change of Chibuzhang Co will be determined and discussed in connection to glacier elevation change later. As shown in Fig. 2b, the horizontal distance between sites A and B is 16 km, and Site A is about 92 km to the central segment of pass 242 over Chibuzhang Co. Figures 2c and d show detailed elevation contours around the two sites. From our visual inspections of optical satellite images on Google Earth, Landsat images, we conclude that site A is close to a glacier terminus and site B is over an icefield. Table 1 also shows the mean glacier slopes along pass 155 around the two sites. The slopes 2.3° (A) and 0.8° (B) are relatively small (< 3°), making it possible to construct a useful time series of glacier elevation changes from altimeters using the dedicated processing procedure given in sections 4.1 and 4.2. Other points along pass 155 have steep slopes and cannot result in glacier changes from T/P altimeters.

ALTIMETER AND SATELLITE IMAGERY DATA
3.1 TOPEX/Poseidon, Jason-1, Jason-2, and Jason-3 Altimeter Data The altimeter data along passes 155 and 242 of the T/P, J1, J2, and J3 altimeters were downloaded from the AVISO website (see section 1). The original elevations from AVISO were improved using the waveform retracking method in this paper (see section 4.1.2). 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-Jason satellites. We used T/P data only from cycles 11 (starting from 1993) to 364 in this paper. The T/P data consist of geophysical data records (GDRs) at 10 Hz (~660 m alongtrack sampling interval) and Sensor Data Records (SDRs) with 64 waveform gates. The nominal tracking gate of a well-positioned tracking window 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 Tanggula Mountains were not within the tracking windows of the J1 radar altimeter. In addition, several problems in the data processing of J1 rendered unreliable data over land (Kuo et al. 2015;Hwang et al. 2016a). In summary, the J1 altimeter data were only 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. 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 retracking gate is 32.5. Finally, the J3 mission was launched in January 2016 and is the fourth of the T/P-Jason satellites. We used J3 data from cycles 1 to 150 in this study. The waveform format of J3 is the same as that of J2.
In this study, 20 Hz SARin observations were used to derive the time series of elevation changes. Off-nadir range correction is applied to the altimeter range (Armitage and Davidson 2014;Abulaitijiang et al. 2015). Since the Cryosat-2 orbit is not exact repeat it is difficult to hit exactly the same spot all the time over rugged terrain like Tanggula Mountains. For all satellites, the height accuracy can be degraded if the radar reflection spot is slightly off the nadir.  In this paper, we will examine height accuracies from Cryosat-2-altimeter observations over Tanggula Mountains and Chibuzhang Co by comparisons with height changes from the T/P-series altimeters. Figure 3 shows the elevations over Tanggula Mountains observed by Cryosat-2 from April 2010 to April 2019, as well as elevations observed by T/P (one cycle) and ICESat. Since the ICESat mission resulted in only one single track over Tanggula Mountains, its elevations were only used for reference in this paper. Note that for lake level changes over Chibuzhang Co, the Cryosat-2 data (GDR-L2) from ESA were used.

SRTM, TanDEM-X, and Landsat-5/7/8 Satellite Data
To assess the altimeter results at sites A and B, we will compare the altimeter-derived glacier elevation changes with those from the following sources: (1) the elevation differences between the DEM from the Shuttle Radar Topography Mission and the DEM from TanDEM-X (Rizzoli et al. 2017), and (2) the area changes from Landsat images. The DEMs from SRTM in the EGM96 orthometric height system are available at grid resolutions ranging from 1" to 15'. In Tibet, the finest grid resolution is 3" in version 2 prior to 2015. TanDEM-X has a finer horizontal resolution than that of the SRTM DEM and the ASTER GDEM (Grohmann 2018). We obtained the needed SRTM DEMs from the USGS SRTM site (https://www.usgs.gov/centers/ eros/science/usgs-eros-archive-digital-elevation-shuttleradar-topography-mission-srtm-1-arc?qt-science_cen-ter_objects=0#qt-science_center_objects). The TanDEM-X DEM was obtained from the Earth Observation Center (EOC) of the German Aerospace Center (DLR) (https:// geoservice.dlr.de/web/dataguide/tdm90/#access). We replaced the geoidal heights in the original TanDEM-X DEM by those based on the EGM96 model to be consistent with the geoid model for SRTM. By using the same geoid model, the geoidal effect is cancelled when computing the elevation difference between these two DEMs over a targeted glacier site. Since our study is concerned with elevation changes, the geoidal model for the altimeter-observed elevations and for satellite-imagery based DEMs must be consistent (EGM96 in this paper). Around sites A and B, we assume that the errors in the SRTM DEM due to co-registration, elevation and radar penetration are a constant and have no effect on the glacier elevation changes determined by differencing the glacier heights from the T/P-Jason altimeters. However, we did consider these errors in the SRTM DEM when determining glacial elevation changes from Cryosat-2 (see section 4.2). Furthermore, we assume TanDEM-X is bias-free over snow and glaciers because TanDEM-X used a short wavelength (X-band) radar.
The Landsat images are from missions Landsat-5, -7, and -8 at 16-day interval. Landsat-5 was launched by NASA on 1 March 1984, and ended on 21 December 2011. Landsat-5 carried the Thematic Mapper to collect images with 7 wave bands ranging from 0.45 to 12.5 µm. Landsat-7/-8 was launched on 15 April 1999 and 11 February 2013, respectively, and continues its operation to date. Landsat-7 provides images with 8 spectral bands and Landsat-8 provides Fig. 3. Elevations (in meter) over part of Tanggula Mountains from the T/P-Jason, Cryosat-2, and ICESat altimeters. A-F are sites for Cryosat-2 determination of glacier elevation change (sites A and B are along T/P pass 155). 11 bands. It is noted that the Scan Line Corrector (SLC) of Enhanced Thematic Mapper Plus (ETM+) onboard Landsat-7 malfunctioned on 31 May 2003, causing strips in the images acquired afterwards. A linear interpolation for the classified image was hence applied in the image processing stage. All Landsat-5/-7/-8 images were downloaded from the Earth Explorer website of the U.S. Geological Survey (USGS; http://earthexplorer.usgs.gov) in the GeoTiff format. In this study, the images covering Chibuzhang Co and Tanggula Mountains are in Landsat-7's frames 139/37 and 138/37 (path/row), respectively, given on the World Reference System-2 grid.

Time Series of Glacial Elevation Change
The flowchart of constructing time series of glacier elevation changes from the repeat altimeter missions is shown in Fig. 4, with details about waveform retracking and quality data selection given in sections 4.1.2 and 4.1.3. First, the orthometric height (glacier elevation) at the footprint of an altimeter radar pulse can be expressed as (Hwang et al. 2016a) Where H g glacier elevation h sat , is ellipsoidal height of the altimeter, R alt is the radar range measurement, N EGM is the geoidal undulation from an Earth Gravitational Model and C is a term containing the following environmental, geophysical and processing corrections: where the first six terms on the right-hand side of Eq.
(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 the T/P satellite and the SGDRs of the Jason series of satellites. The last three processing corrections are related to the radar range correction by waveform retracking (C retrack ; see section 4.1.2), the slope correction (C slope ), and lateral height reduction (C red ), which are explained in detail below. First, the slope correction C slope can be determined as (Fernandez et al. 2015) where R alt is defined in Eq. (1), and a 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 Fig. 2). We found that the C slope values around sites A and B (for all repeat cycles) can easily exceed 1000 m, despite the small slopes around these two sites (see Table 1). Other rugged terrains in Tanggula Mountains can easily have slopes much larger than 2°, which result in unrealistically large C slope values. Because C slope may create great uncertainties in the radar-measured glacier elevations, we decided not to apply this correction in this paper. Instead, we carefully selected 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 thus nearly the same slope corrections (at the same site), thus their slope effects can be reduced or even cancelled when differencing such elevation measurements to determine 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. For the height reduction C red we assume 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 (Lee 2008;Lee et al. 2013). Under this assumption, the height reduction is computed as: where ( , ) DEM 0 0 z m and ( , ) DEM z m are the elevations at the reference point (at latitude 0 z and longitude 0 m , Table 1) and the measurement point (at z, m) from a DEM. The geodetic coordinates (z, m) of measurement points are 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 reduction in Eq. (4). The needed elevations were computed from this SRTM DEM by the bilinear interpolation method.
Over the Tanggula Mountains glaciers, 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-Jason satellites (Davis and Ferguson 2004;Frappart et al. 2006;Hwang et al. 2006;Lee et al. 2013). Frappart et al. (2006) showed 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. (2013) showed that there was no strong correlation between glacier elevation changes and T/P-Jason satellites' Ku-band backscattering coefficients. As such, this study did not consider the effects of elevation change caused by variations in backscattering coefficient.
The waveform retracking for C retrack [Eq.
(2)] and the quality data selection (Fig. 4) are presented in sections 4.1.2 and 4.1.3 for a best possible detail.

Waveform Retracking
The method of waveform retracking for the T/P-Jason altimeters is based on the subwaveform threshold algorithm (Yang et al. 2012). Our subwaveform algorithm first determines a subwaveform of the full waveform that has the maximum correlation with a subset of a Brown waveform (Brown 1977). A Brown waveform can be expressed by the error function with a waveform amplitude, midpoint time and slope of the leading edge and a parameter indicating the decay speed of the trailing edge (Brown 1977;Deng and Featherstone 2006). A reference waveform contains 22 gates and gate values (return powers). Because a full T/P waveform contains 64 gates, a total of 43 correlation coefficients were computed to select the best subwaveform for retracking. Full waveforms of J2 or J3 contains 104 gates, resulting in 83 correlation coefficients, from which the best subwaveform was selected for retracking. Note that our current subwaveform algorithm may be improved using the adaptive approach of Passaro et al. (2014).
A selected subwaveform was used to determine the retracked gate by the threshold algorithm (Davis 1997;Davis and Ferguson 2004). We experimented with three threshold values, 10, 20, and 50%, in the threshold retracking. The use of 50% results in glacier elevation changes that have minimum fluctuations, clear seasonal signals and trends of glacier elevation changes and are well correlated with the trends of glacier area changes (section 4.4). Thus 50% is the adopted threshold value for our final altimeter result presented in this paper. The use of 50% is consistent with the suggestion of Lee (2008). A glacier height correction due to retracking [the term C retrack in Eq. (2)] is the difference between the default gate and the retracked gate, multiplied by a gate-to-meter parameter.
It turns out that only limited waveforms along pass 155 of the T/P-Jason satellites are specular and/or Brownian (Brown-like) and can be retracked to improve the altimeterobserved glacier elevations around sites A and B. As an example, Fig. 5 shows the waveforms from T/P cycle 140 (June 1996; without snow effect to radar) from site A to site B over Tanggula Mountains. The waveforms near sites A and B are specular or Brownian. In particular, waveform No. 10 (specular) is the closest to A and No. 28 (Brownian) is the closest to B. Retracking these two waveforms result in high-quality elevation measurements. Another example is given in Fig. 6 showing J3 waveforms from cycle 121 in May 2019. In Fig. 6, waveforms near site A (No. 20 is the closest) are quasi-specular for effective retracking for improved elevations. However, waveforms near site B (No. 59 is the closest) have multiple ramps, making them not ideal for effective retracking. In Figs. 5 and 6, most waveforms away from sites A and B are very noisy due to off-nadir scattering (or snagging) and rough terrains. Over non open-ocean surfaces, a comprehensive classification of waveforms can be found in Gommenginger et al. (2011). It is also possible to develop an expert system to recover more precise elevations over Tanggula Mountains from the T/P-Jason altimeter measurements.

Quality Data Selection
We selected quality data in three steps as detailed below. Initially we remove erroneous elevation measurements (step 1). In step 2 we select only elevation measurements on the same slope as that at the reference point (Table 1), and in step 3 we fit the selected measurements by a 2nd order surface function while removing outliers. These three steps are explained below.
Step 1: Removing erroneous measurements identifying erroneous elevation measurements by comparing altimeter-measured elevations with those from a SRTM DEM. Following their recommendations, we considered an elevation measurement erroneous if its absolute difference with the 3" SRTM DEM-derived elevation exceeds 150 m.
Step 2: Selecting only elevation measurement on the same slope We selected elevation measurements that are on the same sloping terrain, so that the slope effect can be minimized or diminished [see Eq. (3)]. In addition, glaciers not on the same slope may not have the same rate of elevation change (Forsberg et al. 2017) and they should not be used for a time series that is assumed to have a uniform glacier elevation change rate. As shown in Figs. 2c and d, T/P pass 155 is not parallel to the directions of the steepest descents of the glacier surfaces 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 (Chelton et al. 2001 then altimeter measurement is selected for further processing.
Step 3: Surface fitting and outlier removal The selected altimeter measurements from all repeat cycles in the 1.5-km window (after editing in step 2) are fitted by the following function of time and space (Flament and Remy 2012;Hwang et al. 2016a): where ( , , ) H t i j z m is glacier elevation from repeat cycle j at point i with geodetic latitude (z), longitude (m) [same as those in Eq. (4)] and measurement time (t), V i j is the residual, 0 z and 0 m are the latitude and longitude of the reference point (Table 1), a 0 is the mean elevation, a 1 is the initial rate of elevation change, and a 2a 6 are the coefficients of the function. In Eq. (7), we assume that the altimeter elevation measurements from all repeat cycles must fall into a 2nd order surface that linearly changes with time. The 7 coefficients in Eq. (7) were estimated by the least-squares methods minimizing the weighted sum of the squared residuals. The weight for a raw measurement ( , , ) H t i j z m 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 Eq. (7), for every measurement we examined whether V i j was 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 All the adjusted elevations in a repeat cycle j were then reduced to the elevations at the reference point using the height reduction in Eq. (4): The last step of data processing (Fig. 4) is smoothing. We used the Gaussian filter with a window size of 0.5 years to smooth the time series of glacier elevation changes from the T/P-Jason altimeters. By averaging all such elevations, we computed a representative elevation for repeat cycle j as:

Glacier Elevation Change from Cryosat-2
The method for determining glacial elevation changes from Cryosat-2 is different from the method for the T/P-Jason altimeters because Cryosat-2 does not have an exact repeat period. This method uses a DEM from SRTM for reference elevations and has been widely used in high-mountain Asia (e.g., Kääb et al. 2012;Gardner et al. 2013;Neckel et al. 2014;Chao et al. 2017). The method follows the work by Chao et al. (2017) is given below: (1) The heights from CryoSat-2 (section 3.2) were converted to the orthometric heights relative to the EGM96 geoid, which is the geoid model for SRTM. We removed the outliers in the CryoSat-2-observed glacial elevations using the empirical procedure of Chao et al. (2017).
(2) The SRTM DEM was corrected for the errors due to the universal co-registration (Nuth and Kääb 2011), elevation effect , and radar penetration . The radar penetration effect around the study areas (Fig. 3) is 2.5 m (Liu et al. 2020).
(3) The final CryoSat-2 glacier footprints were determined using the threshold ratio of 2.2, based on Landsat-7 images, Chinese Glacier Inventory (CGI) , and Google Earth images. That is, if the BAND3 and BAND5 indices over a CryoSat-2 footprint satisfies the condition BAND3/BAND5 ≥ 2.2, then this footprint is over glaciers. (4) The width of the data window around sites A, B, C, D, E, F (Fig. 3) is 0.1°. Within the window, the median of the differences between the CryoSat-2 and SRTM elevations is computed. All such medians at the same site form a time series of Cryosat-2-derived glacier elevation changes, fitted by a line that gives an estimated rate of the elevation changes.

Lake Level Change from Altimeter Observations
To discuss the hydrological consequence of Tanggula Mountains' glacier melt in Basins I and III (Fig. 2), we computed lake level changes over Chibuzhang Co from the altimeter data of the T/P-Jason (Pass 242) and Cryosat-2 satellites. Because altimeter waveforms over a calm, open lake surface are mostly specular the accuracies of the altimeter-observed lake levels over Chibuzhang Co will be very high compared with the accuracies over the glacier surfaces of Tanggula Mountains. Our data processing procedure for lake level change is the same as the one used in Hwang et al. (2016b), and is summarized in the following four steps.
Step 1: Selecting altimeter data over Chibuzhang Co First, the lake polygon derived from the GSHHG database (Wessel and Smith 1996) was used to determine Chibuzhang Co's coverage and then a window was set for data selection. According to the experience of Hwang et al. (2016b), 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 used only altimeter data that are 2 km away from Chibuzhang Co's shores.
Step 2: Determining lake elevations from the T/P-Jason altimeters A lake's ellipsoidal height at a given spot of the lake 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 (section 4.1.2) to determine the range corrections for the lake elevation measurements. However, a 20% threshold value was used in the retracking over Chibuzhang Co, instead of 50% that was used over glaciers. The use of 20% is based on the result of Hwang et al. (2016b) 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 no further outliers were 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 segment of pass 242 over Chibuzhang Co is short (Fig. 2b) and the selected altimeter data are near the lake center, we neglected the geoidal height differences over this lake segment of pass 242 when averaging the elevation measurements to obtain the mean elevation.

Glacier Area Change from Landsat Imagery
For comparison with the altimeter-derived glacier level changes, we determined glacier area changes around sites A and B using the Landsat images described in section 3.3. In this paper, the annual area change is defined as the smallest detected area change each year in the glacier-covered areas near site A or B. The steps for determining glacier area changes are described below (Chander et al. 2009); see also the sample images and plots in Figs. 7a -c.
Step 1: Computing the Top-of-Atmosphere (ToA) reflectance We selected cloud-free images over Tanggula Mountains to compute pixel-wise glacier coverages (Fig. 7a). For each pixel, first we computed the Top-of-Atmosphere (ToA) reflectance using (Cea et al. 2007) where L m is the radiance, gain and bias are the gain and bias for the Green or the mid-infrared band, and DN is the digital number of the pixel. Then we determined where t m is the ToA reflectance (unitless), d is the distance between the Earth and Sun, ESUN m is the mean solar exoatmospheric irradiance, and s i 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, Fig. 7b) using where G t and M t are the reflectances in the green and midinfrared bands from Eq. (12). Following the recommendation of Cea et al. (2007), we decided that if NDSI ≥ 0.4, then the pixel is over glacier. The resulting image was compared with the original image to ensure the glacier coverage is correct. As shown in Fig. 7c, for each image the glacier edge is detected using NDSI, and the edge defines the glacier cov-erage at the image-acquisition time. We selected the smallest glacier-coverage area during the dry season of a year to construct a time series of glacier area changes that is used to assess the glacier level changes from the altimeters (see section 5.2). Figures 8a and b show the altimeter-derived glacier elevation changes at sites A and B, both suggesting steady glacier elevation declines. Compared with previous studies such as Lee et al. (2013), we used a stricter selection criterion (see Fig. 3) resulting in less elevation observations from the T/P-Jason altimeters. In addition, there are only few qualified elevation measurements from J2 points at sites A and B. Table 2 shows the rates of glacier elevation change and the rates of lake level change (the discussion about lake level change will be presented in section 6). The rates of glacier elevation change from T/P at sites A and B are close at about -3 m year -1 over 1993-2002. Table 2 shows that only 22-27% of the raw T/P altimeter measurements were selected (or qualified) at the two sites. The percentages of the usable J3 data are slightly lower than those  Table 2. Altimeter-derived rates of glacier and lake elevation changes.

Note: a: For glacier elevation changes at A and B, this is amplitude of annual oscillation from T/P measurements.
of T/P. The low percentages (about 1/4 of the raw data) of the usable T/P and J3 elevation measurements are mainly due to the quality data selection procedure given in section 4.1.3 (see Fig. 4). In contrast to the low percentages (25%) of the usable T/P data over the two glacier sites, the usable percentage of lake level measurements over Chibuzhang Co is 69% (Table 2). Figures 8a and b also show the glacier elevation changes from Cryosat-2 at sites A and B, which were consistent with the elevation changes from J2 and J3. Figure 9 shows the time series of glacier elevation changes at sites A -F with and without the corrections for the SRTM DEM (section 4.2). The corrections result in smaller oscillations in elevation changes from Cryosat-2. Figure 9 suggests that during 2010-2019 the largest declining rate of glacier elevations occurred at site A, followed by site B. Over sites D -F, the glacier declining rates are smaller and perhaps this is the result of the altitude effect: the higher the glacier, the less it melts (sites D -F are above 6000 m).
Despite the small numbers of qualified elevation measurements from J2 and J3 over 1993-2020 at site A, we estimate the rate of glacier elevation changes from J2 and J3 here. The resulting rate over 1993-2020 is -1.87 ± 0.05 m year -1 , compared with -1.12 ± 0.050 m year -1 from Cryosat-2 ( Fig. 9). At site B, the rate from T/P and J3 over 1993-2020 is -1.46 ± 0.05 m year -1 (no J2 data are usable for this rate), compared with the rate of -0.65 ± 0.34 m year -1 from Cryosat-2. The rates given in Table 2 and Fig. 9 shows that the glacier thinning rates at sites A and B have decelerated from -3 m year -1 over the period of 1993-2002 (from the T/P measurements only) to around -1 m year -1 over the period of 2002-2020 (from the J2, J3, and Cryosat-2 measurements). Despite some inconsistencies in these results, the rates of melting are in general the same. The rates vary over space and time, and over different periods. In addition, two notable features of the time series in Figs. 8a and b are: (1) there was no winter glacier elevation high in 1997-1998 at site A, and (2) there were large glacier elevation fluctuations at site B in the spring of 1998 (from T/P) and in the summer and fall of 2017 (from J3). A possible link of such abnormal glacier elevation changes to the 1997-1998 El Niño and to lake changes of Chibuzhang Co will be discussed in section 6. Figure 10 shows a satellite image around site A, indicating that this site is over a deep valley. The satellite images from Google Earth over 1984-2016 around site A show that glaciers near the footprints of T/P at site A nearly vanished by 2016. That is why the J3-observed glacier elevation changes in Fig. 8a were relatively flat: the observed elevations near 2016 are almost over a non-glacier surface, thus showing no continuous elevation declines, which should occur over a surface with melting glacier. By contrast, site B is over a large icefield with fluctuating elevations affected by glacier melt. Again, at site B, there were large elevation fluctuations in 1997-1998 and 2017 (Fig. 8b), which require more investigations about the potential causes. One likely cause of the large fluctuations in glacier elevation at site B is the rising temperature and increased precipitation that may have roughened the glacier surfaces here. Such roughened surfaces may have led to contaminated waveforms that cannot be repaired by the subwaveform retracking (section 1.1.2). Large elevation oscillations were also found in the altimeter-derived elevation time series over Bering Glacier in Alaska (Lee et al. 2013).

Comparison with DEM Differences and Landsat-Derived Area Changes
Our first assessment of the altimeter result (section 5.1) is comparison with the direct elevation differences between the SRTM and TanDEM-X DEMs (section 3.3). Figure 11 shows the elevations from these two DEMs and the elevation differences around sites A and B. Figure (Rizzoli et al. 2017), which range from December 2010 to early 2015 (note that the image acquisition time for the SRTM DEM is from 11 to 22 February 2000; no radar penetration effect is applied here). At site A, the DEM-derived rate of about -3 m year -1 (based on the time span from 2000 to 2011) is larger than (in absolute magnitude) the rate of -1.87 m year -1 derived from the T/P, J2, and J3 altimeters (1993-2020; Table 2). This consistency in rates also happens at site B: the DEM-derived rate of about -1 m year -1 is close to the rate of -1.46 m year -1 from the T/P and J3 altimeters (1993-2020, Table 2). This assessment shows that the altimeter-derived rates of glacier elevation change at sites A and B are in general consistent with those from the direct elevation differences between the TanDEM-X and SRTM DEMs.
The second assessment, which is an indirect assessment, is based on the Landsat-derived area changes (section 4.4). Figures 12a and b show the time series of annual minimum glacier area around sites A and B. Table 3 shows the glacier area changes over different periods from 1986 to 2019 from the Landsat imagery. The glacier area changes near sites A and B for over 1986-2019 in Table 3 were negative and the total changes are around -2.20 km 2 in this period. In the processing of Landsat-series images, we calculated the minimum of the annual snow/ice area from as many cloud-free images as possible. We assume that snow cover did not affect this minimum area. However, some Landsat images failed to deliver results over Site A due to cloud covers that made it difficult to identify snow coverages using Eq. (13).
At site A, the glacier areas declined almost monotonically, except the slight increases after 2016. At site B, the pattern of glacier areas changes exhibits unstable oscillations between 1996 and 2010 and the mean rate turned positive after 2016. However, the mean rate of glacier area changes at site B over 1986-2019 was negative. In addition, the glacier areas at site B fluctuated more rapidly than those at site A. The Landsat imagery delivers only yearly minimum glacier areas without seasonal area changes. In contrast, the T/P-Jason 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 Figs. 8 and 12, we can draw the following conclusions: (1) the glacier elevation changes at site A from the T/P-Jason altimeters were positively correlated with the Landsat-derived glacier area changes; (2) at site B, the glacier elevation changes from T/P over 1993-2002 were positively correlated with the glacier area changes; 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 (Fig. 8b).
In summary, the second assessment using the Landsatderived area changes has largely confirmed the glacier thinning at sites A and B given in the altimeter result. However, it is the first assessment (direct comparison with the SRTM-TanDEM-X DEM difference) that provides the direct validation of the altimeter-detected glacier thinning sites A and B.

GLACIER LEVEL CHANGE AND LAKE LEVEL CHANGE
In this section we discuss how the glacier losses over Tanggula Mountains may affect lake level changes using elevation changes observed by the altimeters (sections 3.1 and 3.2). Chibuzhang Co and its sister lake Dorsoidong Co are located in the west of Tanggula Mountains. The two lakes are glacier-fed lakes and receive glacier melts from Tanggula Mountains (Qiao 2010;Song and Sheng 2016;Chao et al. 2017). Note that most of the time the middle Lake in Fig. 2a has the same level as Chibuzhang Co and is part of the latter. In this paper, these two lakes are named as one lake-Chibuzhang Co (Figs. 2a and b). Figure 13 shows the lake level changes over Chibuzhang Co and the glacier level changes from the altimeters. The lake level changes from the T/P-Jason altimeters and Cryosat-2 are consistent over 2010-2020. Both lake level and glacier level changes in Fig. 13 show clear seasonal oscillations, but they have the opposite trends. In general, the glacier elevation of Tanggula Mountains 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. Lake level changes can be caused by many factors, such as glacier melting, precipitation, evaporation, human activity, and tectonic activity. The lake level highs of Chibuzhang Co in summers are partially caused by glacier (e) (f) Fig. 9. Glacier elevation changes at sites A to F from Cryosat-2.  Table 3. Landsat-derived glacier area changes (km 2 ) over different periods. meltwater from Tanggula Mountains and partially by rain in Basin I (Fig. 2a). As shown in Fig. 2b, the distance between T/P pass 242 (lake center) and the nearest equilibrium line of 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 (Fig. 2b), rather than Basin I, site A's short distance to Basin I implies that the glacier elevation changes observed at site A (Basin III) should be highly correlated with those in Basin I. That is, a T/P-Jason-observed glacier melt at site A might signify a glacier melt in Basin I that in turn can lead to the lake level rise seen in Fig. 13, particularly after 2003.
The result in Fig. 13 shows that the rate of lake level rise in Chibuzhang Co is 0.06 ± 0.01 m year -1 over 1993-2002 (T/P only), increasing to 0.18 ± 0.01 m year -1 over 2003-2020 (all altimeters, Table 2). That is, the rate in the latter period (2003-2020) is almost three times the rate in the former period (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002), implying an accelerated supply of source water to Chibuzhang Co. This result supports the results of Tseng et al. (2016) and Jiang et al. (2017) who proposed that global warming has increased the meltwater from Tanggula Mountains, which then led to the connection between Chibuzhang Co and Dorsoidong Co sometime between 2006 and 2008. Because the mean elevation of the original Chibuzhang Co was higher than that of Dorsoidong Co before 2006, it is suggested that the persistent meltwater from Tanggula Mountains has steadily increased the level of Chibuzhang Co to flood Dorsoidong Co. Eventually the two lakes were connected.
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, Figs. 2a and b) and Chibuzhang Co, respectively. This implies that the middle lake is about 23% of Chibuzhang Co in area. The channel from the middle 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 during the mission time span of T/P (1993-2002) seen in Fig. 13.
It is likely that increased precipitation and temperature-induced glacier melt have contributed to the rising level of Chibuzhang Co. This issue has been discussed in several studies. For example, Wu et al. (2015) and Chao et al. (2017) showed that precipitation increased at the rates of 2.714 and 2.285 mm year -1 over 1988-1999 and 1966-2013, respectively, at the Tuotuo River meteorological station near Tanggula Mountains. In addition, Chao et al. (2017) showed an increasing rate of 0.032 °C year -1 in temperature over 1979-2007 at the Tuotuo River station.
The mass balance of Xiao Dongkemadi glacier (near Tanggula Mountains) has been studied by Wu et al. (2015), who found that the elevation of Xiao Dongkemadi's snow equilibrium line increased in 1997, but decreased suddenly in 1998. Pu et al. (2008) speculated that these dramatic elevation oscillations in the snow equilibrium line were caused by the 1997-1998 El Niño. The lake-level time series in Fig. 13 contains anomalously high values in the summer of 1998. This high may be associated with the sudden elevation rise in the snow equilibrium line in the summer of 1998, which implies more meltwater to enter Chibuzhang Co. In addition, Figs. 8a and b show several anomalous glacier elevation changes around 1997-1998. First, at site A (Fig. 8a) 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 elevation should be lower). Furthermore, at site B (Fig. 8b) there was a large glacier elevation drop in the spring of 1998, which was coincident with the sudden increase in the lake level of Chibuzhang Co (Fig. 13). The difference between the patterns of glacier elevation change at sites A and B around 1997-1998 could be attributed 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 the noise levels of altimeter radar ranging.

CONCLUSION
This paper shows that it is possible to measure seasonal, inter-annual and secular elevation changes over inaccessible, small-patched mountain glaciers by the T/P-Jason and Cryosat-2 altimeters. Cyrosat-2 detected glacier elevation changes consistent with those from J2 and J3, despite the fact that Cyrosat-2 is not in a repeat orbit. The Cryosat-2-observed glacier elevations show the altitude effect of glacier melting: the higher the glacier, the less it melts. The altimeter-observed glacier thinning in Tanggula Mountains is confirmed by the direct elevation differences between the DEMs from the missions TanDEM-X and SRTM, and indirectly by the glacier area losses from Landsat imagery.
The glacier thinning was correlated with the rising lake levels around Tanggula Mountains. In particular, the altimeter-observed, accelerated rate of lake level rise in recent years may signify growing precipitation and more glacier meltwater in this region. The rising lake levels have led to a new lake formed by Chibuzhang Co and Dorsoidong Co. The new lake could have modified the original wetland system and the ecological system. In contrast to the accelerated lake level rise of Chibuzhang Co, the thinning rates of the glaciers around sites A and B have decelerated from the period of 1993-2002 to the period of 2003-2020, likely due to the persistent glacier losses here that no longer can sustain large glacier thinning rates in recent years.
This study shows the potential of a nadir-looking altimeter in providing an altimeter-based, virtual glacier station. With the algorithm used in this paper, only about 1/4 of the raw T/P elevation measurements at s sites A and B can be used for a virtual station, and 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. Thus, increasing the volume of quality altimeter data and improving the radar ranging accuracy over mountain glaciers remain two challenging issues when implementing the concept of altimeter-based virtual glacier station.