**3. Methods**

Since the ICESat-1&2 data are laser point data and the spatial distribution of the data points is not uniform in space, the ICESat-1&2 data therefore cannot provide pixel-wised information covering the entire surface of glacier like the ZY-3 stereo image pair data. The spatial distribution of ICESat-1&2 data points on each glacier differs due to different satellite overpassing tracks and different sensor laser beam trajectories. In addition, due to the different atmospheric conditions and different radiation energy received by glaciers in different elevation and orientations, the glacier mass balance will differ in different aspects at the same altitude. Based on the above reasons, this paper proposed an analysis method, namely "elevation-aspect bin analysis method", to investigate glacier elevation changes by considering the glacier terrain elevation, slope aspect and the corresponding glacier areas.

Our research workflow consists of: (1) extraction of the elevation differences between the ICESat-1&2 data and the SRTM DEM data in the ICESat-1&2 footprint scale; (2) extraction of glacier elevation difference in 1◦ × 1◦ grids; (3) estimation of the change of glacier elevation in sub-regions or river basins; and (4) uncertainty analysis for the quantitative assessment of the influences of data errors, processing errors and other factors on the

results. Figure 3 shows the specific workflow, and the method details are described in the following sections.

**Figure 3.** Flowchart of glacier elevation change extraction and analysis using ICESat-1&2 data and SRTM DEM data.

### *3.1. Extraction of Elevation Difference in the ICESat-1&2 Data Footprints*

Due to the differences in reference ellipsoid for projection and differences in vertical reference system between ICESat-1&2 data and SRTM DEM data, we first need to convert the elevation datum to extract elevation change information and eliminate outliers. In this study, bilinear interpolation method was used to extract the elevation information of SRTM DEM data in the footprint of the ICESat-1&2 data. Two corrections were applied before the difference between the elevation from ICESat-1&2 data and from SRTM DEM data can be calculated: (1) The Topex/Poseidon ellipsoid used by ICESat-1 data was converted to the World Geodetic System 1984 (WGS84) ellipsoid used by ICESat-2 data and SRTM DEM data [49,50]. (2) The Earth Gravity Model (EGM) 1996 used by the SRTM DEM geoid elevation was converted to EGM2008 used by ICESat-1&2 data [35].

The glacier elevation change in each year, Δ *Hyr* (in meter), at each footprint point can then be calculated by:

$$
\Delta H\_{yr} = H\_{ICESat, yr} - H\_{SRTM} + P\_{SRTM} \tag{1}
$$

where *HICESat*,*yr* (m) is the elevation from ICESt-1&2 data in a year "*yr*" after ellipsoid correction; *HSRTM* (m) is the elevation from SRTM DEM data in 2000 after leveling correction and bilinear interpolation within the footprint of ICESat-1&2 data point; *PSRTM* (m) is the penetration depth of the SRTM DEM data in the glacier region, we set *PSRTM* as 2.4 m in this study according to [51].

To reduce the influence of slope, ICESat-1&2 data points in slopes greater than 30◦ were excluded according to [35,52]. To remove outliers affected by detector saturation, only ICESat-1 GLAS14 data flagged with class "0" or "1" in the data quality layer "satCorrFlg" were selected according to [53]. For ICESat-2 data, we only selected the data flagged with category "0" in the data quality layer "atl06\_quality\_summary" to ensure the data quality according to [48]. Additionally, we excluded the data points with elevation differences greater than 100 m between ICESat-1&2 data and the SRTM DEM data to eliminate the influence of outliers that were possibly introduced by cloud interference [35].

#### *3.2. Extraction of Glacier Elevation Difference in 1*◦ × *1*◦ *Grids*

To reduce the uncertainty caused by the uneven spatial distribution of the ICESat-1&2 data points in each elevation and slope aspect on the results, the core idea of the "elevationaspect bin analysis method" proposed in this paper is to calculate the glacier elevation change in each 1◦ × 1◦ grid according to the distribution of ICESat-1&2 data points in different elevations and slope aspects and respective glacier areas. The whole HMA is divided into 1◦ × 1◦ grids. In each 1◦ × 1◦ grid, the elevation is divided into elevation bins (denoted by *i*), each elevation bin is divided into 8 aspect bins (denoted by *j*) with an interval of 45◦.

The elevation difference for each elevation bin in a 1◦ × 1◦ grid is calculated by:

$$
\Delta H\_{yr,\,grid}(i) = \sum\_{j=1}^{8} \left( F\_{yr,\,grid}(i,j) \cdot \Delta H\_{yr,\,grid}(i,j) \right), \ (j = 1, 2, \dots, 8, \text{number of aspect bins}) \tag{2}
$$

where, *Fyr*,*grid*(*<sup>i</sup>*, *j*) and <sup>Δ</sup>*Hyr*, *grid*(*<sup>i</sup>*, *j*) are the fraction of glacier area and the median of glacier elevation difference (between a target year and 2000) in aspect bin *j* of the elevation bin *i* in a 1◦ × 1◦ grid.

The elevation difference between a target year and the reference year 2000 in a 1◦ × 1◦ grid, <sup>Δ</sup>*Hyr*,*grid*, was the area-weighted average of elevation difference for all elevation bins and calculated as:

$$\Delta H\_{yr,grid} = \sum\_{i=1}^{N} \left( F\_{yr,grid}(i) \cdot \Delta H\_{yr,grid}(i) \right), \ (i = 1, 2, \dots, N, \text{ number of elevation bins}) \tag{3}$$

where *Fyr*,*grid*(*i*) is the fraction of glacier area of the elevation bin *i* in a 1◦ × 1◦ grid. Following the method in [35], we repeated the calculation using Equations (2) and (3) by setting 6 different intervals of elevation bin, i.e., 200 m, 300 m, 400 m, 500 m, 600 m and 700 m, and used the average of the six results as the final glacier elevation difference in each 1◦ × 1◦ grid. ICESat-1 and ICESat-2 data from September to November of each year in 2003–2008 and 2018–2020, respectively, were used for annual change calculation.

Monthly glacier elevation difference between all months in 2019–2020 and the reference year 2000 in each 1◦ × 1◦ grid was calculated in the same way using the "elevation-aspect bin" areaweighting method.

#### *3.3. Estimation of Annual and Monthly Change of Glacier Elevation in Sub-Regions*

The annual glacier elevation changes in a sub-region or a river basin, <sup>Δ</sup>*Hyr*, *subregion*, was calculated by the area-weighted average of glacier elevation differences in all 1◦ × 1◦ grids in a sub-region or river basin as:

$$
\Delta H\_{\text{yr, subreglou}} = \sum\_{k=1}^{M} \left( F\_{\text{yr}}(k) \cdot \Delta H\_{\text{yr, gydd}}(k) \right), \text{ ( $k = 1, 2, \dots, M$ , number of grids in a region)} \tag{4}
$$

where *Fyr*(*k*) is the ratio of glacier area of each 1◦ × 1◦ grid over the total glacier area in either a sub-region or a river basin.

Following the method used in Kääb et al. [54], the rate of change of glacier elevation over multiple years was calculated by applying robust linear regression to the time series of <sup>Δ</sup>*Hyr*, *subregion*.This way, the rate of change in glacier elevation between 2003–2008 and between 2003–2020 were obtained using ICESat-1 and ICESat-2 data, respectively.

### *3.4. Error Analysis*

It is difficult to directly assess the uncertainty of ICESat-1&2 and SRTM DEM data using ground data due to insufficient in situ measurements of glacier elevations. Factors affecting the accuracy of the results include the vertical height error of SRTM DEM data, the penetration depth error of SRTM DRM data in the glacier area and the error of ICESat-1&2 data itself. In this study, the uncertainty assessment method of glacier monitoring by Wang et al. [35] was adopted, as shown below:

$$
\sigma\_{DH} = \sqrt{\sigma\_{std}^2 + \sigma\_{dh}^2} \tag{5}
$$

where *σDH* is the glacier elevation change error (m), *σstd* is the standard deviation of the result in different elevation bins (m). The *σdh* includes the SRTM DEM vertical deviation (<16 m), ICESat-1&2 error (centimeter and decimeter scale) and radar penetration error (meter scale). In this study, the absolute error of each ICSat-1/2 footprint point is set to 20 m.

The uncertainty of glacier elevation change rate is obtained by

$$
\sigma\_{\rm DHI/dt} = \sqrt{\sigma\_{\rm flt}^2 + \sigma\_{\rm spat}^2 + \sigma\_{\rm temp}^2 + \sigma\_{\rm hlas}^2} \tag{6}
$$

where *<sup>σ</sup>DH*/*dt* is the glacier elevation change rate error, *<sup>σ</sup>fit* is the linear fitting error, *<sup>σ</sup>spa<sup>t</sup>* and *<sup>σ</sup>temp* are spatial and temporal sampling errors, *σbias* is the comparison bias and unknown system uncertainty of crustal uplift [55]. We followed Wang et al. [22] and set *<sup>σ</sup>spat*, *<sup>σ</sup>temp* and *σbias* to 0.06 m/year in this study.
