Contribution of ground ice melting to the expansion of Selin Co (lake) on the Tibetan Plateau

. Selin Co, located within permafrost regions surrounded by glaciers, has exhibited the greatest increase in water storage among all the lakes on the Tibetan Plateau over the last 50 years. Most of the 15 increased lake water volume has been attributed to increased precipitation and the accelerated melting of glacier ice, but these processes are still not sufficient to achieveclose the water balancebudget with the expansion of Selin Co. Ground ice meltwater released by thawing permafrost due to continuous climate warming over the past several decades was regarded as another source of lake expansion. This study presented the first attempt to quantify the water contribution of ground ice melting to the expansion of Selin 20 Co by evaluating the ground surface deformation. We monitored the spatial distribution of surface deformation in the Selin Co basin using the SBAS-InSAR technique and compared the results with the findings of field surveys. Then, the ground ice meltwater volume in the watershed was calculated based on the cumulated settlement. Finally, this volume was compared with the lake volume change during the same period, and the contribution ratio was derived. SBAS-InSAR monitoring during 2017–2020 illustrated 25 widespread and large subsidence in the upstream section of the Zhajiazangbu subbasin, where widespread continuous permafrost is present. The terrain subsidence rate was normally between 5 and 20 mm/a, indicating rapid ground ice loss in the region. The ground ice melted water reachedwas released at a rate of ~57.4×10 6 m 3 /a, and the rate of increase in lake water storage was ~485×10 6 m 3 /a during the same period, with ground ice meltwater contributing 11.8~12% of the lake volume increase. This study is especially 30 helpful forcontributes to explaining the rapid expansion of Selin Co and equilibrating the water balance at the watershed scale. More importantly, the proposed method can be easily extended to other watersheds underlain by permafrost and to help understand the hydrological changes in these watersheds.


Introduction 35
More than 1200 lakes on the Tibetan Plateau (TP) span an area exceeding 1 km 2 , and the total lake area is greater than 46000 km 2 . The water in most of these lakes is more or less connected with widely distributed glaciers and permafrost. Recent studies have indicated that most of the lakes on the TP have manifested extensive changes . In particular, Selin Co (also known as Siling Co, Serlin Co, and Serling CoCo is lake in Tibetan) exhibited the greatest increases in both lake 40 area and water storage: its lake area expanded by ~40% from ~1700 km 2 in 1972 to ~2400 km 2 in 2020, and its water storage increased by 80% from 309.4×10 8 m 3 in 1972 to 558.4×10 8 m 3 in 2017 . Its lake area surpassed that of Nam Co in the early 2000s Bian et al., 2010); consequently, and Selin Co is now the second largest saltwater lake in China. Such rapid changes in Selin Co have significantly affected the regional environment and have thus attracted substantial interest 45 within the scientific community.
The entire Selin Co watershed covers a drainage area of 4.4×10 4 km 2 , 18 times the lake surface. Accordingly, a lake water level increase of by ~0.2 m/a corresponds to at least 1 cm/a of water collected uniformly over the whole watershed area (neglecting evaporation). The entire watershed hosts 299 glaciers with a total area 50 of 369.7 km 2 and ice reserves of 27.9 km 3 based on the second Chinese glacier inventory ; additionally, according to the new estimation by and an ice volume of 21.8 km 3 , the glacier volume reaches 21.8 km 3 in the watershed. .  mapped the permafrost distribution on the Qinghai-Tibet Plateau (QTP) based on freezing and thawing indices from Moderate Resolution Imaging Spectroradiometer (MODIS) land surface temperatures and 55 validated their map using various ground-based datasets; according to this permafrost map, the permafrost area covers ~1.3×10 4 km 2 , accounting for 30.2% of the watershed. Continuous permafrost and seasonally frozen ground exist, with the former being widespread mainly in the northern part of the watershed .   estimated the ground ice storage on the QTP based on the ice content distribution characteristics from 164 drill core records deeper than 15 m, the above 60 permafrost distribution map, a map of the Quaternary sedimentary types, and a permafrost thickness map; based on their map, the ground ice volume in the watershed reaches 132.3 km 3 , approximately five times the glacier ice volume in the Selin Co watershed.
To reveal the potential reasons for the abovementioned rapid expansion of lakes on the TP, an accurate 65 estimation of the basin water balance is urgently needed . Glacial meltwater, the thawing of permafrost, precipitation (including snow) and changes in evapotranspiration all contribute to lake recharge. A model simulation of endorheic basins on the TP showed that increased net precipitation contributed the majority of the water supply (~70%) to the increased lake volume . In addition, recent research has revealed that glacial meltwater has 70 contributed ~ 10% of the total water input to Selin Co since the 1970s Tong et al., 2016;. The weakening of lake evaporation during 1972-2010 due to decreasing wind speeds also contributed to the accelerated expansion of Selin Co to some extent, but this contribution was reported to be very small .

75
Significant permafrost degradation has been observed on the TP under the impacts of the warming climate.
The monitoring of ten boreholes on the TP revealed that from 1981 to 2018, the active layer thickened at an average rate of 19.5 cm per decade; moreover, this thickening trend has been accelerating in recent years . In the meantime, different permafrost regions across the TP experienced thaw settlements Chen et al., 2022). The ice content within the uppermost layer of permafrost is typically 80 higher than the saturated water content after this permafrost layer thaws; hence, the thawing of this layer might result in terrain settlement . The terrain settlement was attributed to the melting of ground ice from the ice-rich permafrost layer just below the permafrost table and the further release of this water into the hydrological cycle . 85 The melting of ground ice released water into the hydrological cycle Jin et al., 2022), hence permafrost degradation might be a possible source of the water that is causing lakes to expand , but this contribution has yet to be quantified. By separating hydrographs using isotopes, some researchers calculated the contributions of meltwater from ground ice to surface water runoff and found that 90 the contribution ranged from 13.2% to 16.7% in the source region of the Yellow River , and that to nearby thermokarst lakes in the Beiluhe region reached ~61.3% . The amount of meltwater from permafrost degradation was also modeled by multiplying the active layer thickening rate and the average ground ice content .

95
Interferometric synthetic aperture radar (InSAR) analysis can exploit the phase changes in SAR signals to determine relative surface displacements on the order of millimeters to centimeters. Accordingly, InSAR monitoring can detect the terrain subsidence triggered by the thawing of ice-rich permafrost, e.g., Liu et al., 2012a). (e.g., Liu et al., 2012a). Permafrost regions 100 with higher ground ice contents have been shown to produce greater terrain subsidence on the TP . The subsidence rate is less than 5 mm/a in the northwestern TP, where the climate is cold and dry , while it is up to 30 mm/a on the ice-rich Eboling Mountain in the northeastern region of the TP .

105
In this study, we conducted the first attempt to quantify the contribution of ground ice melting to the expansion of the Selin Co watershed through surface settlement. It is well known that the permafrost layer just below the permafrost table always contains ground ice higher than 50% in volume Mackay, 1983;. Therefore, we assumed that the amount of surface settlement would release the same amount of ground ice caused by compressing the thawing ice-rich 110 permafrost layer. The spatial distribution of surface deformation was derived through InSAR analysis and compared with the findings of field surveys. Finally, the ground ice meltwater volume in the watershed was compared with the Selin Co volume increase, and the contribution ratio was determined.

Study area 115
The Selin Co watershed, which is in the transitional zone between the Indian monsoon and the westerlies over the TP, is characterized by a cold semiarid monsoon climate with a mean annual temperature of approximately 0°C and average precipitation of ~350 mm (Tong et al., 2016). From 1979 to 2017, the average annual temperature increased at a rate of 0.049°C/a, and the average annual precipitation increased at a rate of 4.65 mm/a, with the main increase occurring after the mid-1990s . 120 This watershed basin extends over a total area of approximately 4.4×10 4 km 2 , and the rivers and lakes within the basin are connected, forming an inland lake group. Ngoin Co is located south of Selin Co, while Wuru Co & Qiagui Co are located west of Selin Co; their locations are marked as ② and ③, respectively, in Fig.   1. The main rivers draining into the lake are the Zhajiazangbu (from the north), the Zhagenzangbu and 125 Alizangbu (from the west), and the Boquzangbu (from the east) (Tong et al., 2016), as shown in Fig. 1. The details of the four major subbasins and the inflows of the runoff from these rivers into the lake are listed in Table 1. The Zhajiazangbu river, with a length of 409 km, originates from the Tanggula Mountains and enters Selin Co from the north. Widespread continuous permafrost occurs mainly in the northern part of the basin, whereas sporadic permafrost and seasonally frozen ground are found in the central and southern parts of the 130 watershed. Among the four subbasins, the Zhajiazangbu has the largest permafrost distribution (10667 km 2 or 66.2% of the subbasin area), followed by the Zhagenzangbu (1967 km 2 or 12.3% of the subbasin area) .

Field data 145
We conducted a field investigation of the permafrost in the study area in October 2019. Seven boreholes deeper than 20 m were drilled, and ground-penetrating radar (GPR) surveys were carried out in the upstream section of the Zhajiazangbu subbasin (locations are marked in Fig. 1). The survey was performed at the end of the thawing period, allowing us to estimate the maximum thawing depth and the location of the permafrost table. The permafrost table was estimated from the cores of the boreholes, and the development of ground 150 ice was described as well. Descriptions of the boreholes are provided in Table 2, and field photographs of the cores at borehole sites SLC01 and SLC04 are shown in Fig. 2.
Twenty 50-m-long GPR profiles were carried out within or around the catchment, fourteen of which were within the basin. Based on the fourteen GPR profiles carried out within the basin, the permafrost table was 155 at a maximum depth of 4.1 m, a minimum depth of 2.0 m, and an average depth of 3.2 m. The volumetric soil water content in the active layer ranged from a maximum of 46.4% to a minimum of 14.7%, with an average of 22.6% based on our interpretation of the GPR data.

Sentinel-1 SAR images
Sentinel-1 (S1) C-band SAR images (https://scihub.copernicus.eu/) were used to monitor the surface deformation. S1, a constellation comprising two satellites launched in 2014 (S1A) and 2016 (S1B) that orbit the Earth, is a C-band SAR mission that was developed and is operated by the European Space Agency (ESA) 170 as part of the Copernicus program. VV-polarized Level 1 single look complex (SLC) images acquired in interferometric wide-swath (IW) mode in descending orbit were used in the study. The study area is covered by two orbits of S1 acquisitions, the details of which are shownlisted in Appendix Table 3A1. In total, 95 acquisition dates for descending orbit 48 and 100 acquisition dates for descending orbit 150 from September 2017 to December 2020 were processed. 175 Table 3 Information on the S1 images processed to monitor the surface deformation.

Orbit 48
Orbit 150 was used to calculate the slope, and remove the topographic phase and implement geocoding during InSAR processing.
Cloud-free images acquired during autumn or summer were selected, and the acquisition date and path/row of each image are listed in Appendix Table 4A2. 190 "p" denotes the path, "r" denotes the row of the frame, and "D" denotes the acquisition date.

Workflow
An overview of the methodology employed in this study is illustrated in Fig. 3. The main steps are 195 summarized as i) retrieval of lake water storage changes, ii) retrieval of deformation time series, iii) estimation of the ground ice meltwater volume from the long-term deformation rate, and iv) calculation of the ratio of the water volume contributed by permafrost ground ice melting to the increase in lake water storage. The processing steps are described in detail in the following subsections.

Lake water storage changes 205
First, we retrieved water surface elevations through the ICESat-2 ATL13 laser altimetry product by extracting data values within each lake's area from 13 October 2018 to the end of 2020. For each year, the mean annual elevation was calculated from all data acquisitions during the year. The lake water level change rate was then estimated by a linear trend model.

210
Then, we extracted the lake area extent using Landsat 8 OLI images in 2018, 2019, and 2020. Three bands (NIR, SWIR1, and SWIR2) were stacked to extract the lake area because lake water has extremely low reflectance in these bands and thus is easily separated from other land cover types. K-means clustering was then applied to classify the land cover types in the stacked image into two classes: water and other. After that, a 3×3 window size majority filtering was applied as a postprocessing step, and we then converted the 215 classified image into a shapefile and eliminated the spots outside the inlet and outlets of the lakes. The lake areas in 2015, 2016 and 2017 were obtained from a previous study (Zhang, 2019); we also modified the lakes to have the same inlets and outlets. For lakes with irregular areas, their storage was approximately calculated according to the volume of a circular platform; thus, the change in the volume of the lake was calculated from the difference between the volumes of two circular platforms , as described in Eq. 220 (1): where ΔV is the change in lake water storage, 1 and 2 represent the lake areas during two periods (e.g., the lake areas in 2018 and 2020), and ( 2 − 1 ) represents the water surface elevation change between the two periods. 225

SBAS-InSAR processing
In this study, we implemented the small baseline subset (SBAS) time series InSAR analysis (SBAS-InSAR) technique, which deploys multiple master datasets to minimize the effects of spatiotemporal decorrelation Lanari et al., 2004;Usai, 2003) by selecting interferograms with small spatial and 230 temporal baselines. Thus, this technique is suitable for permafrost environments prone to strong spatiotemporal decorrelation.
The Selin Co basin is covered by two orbits of S1 images. We processed the SBAS-InSAR analysis individually for each orbit. SBAS-InSAR processing in this study contains three main steps: i) InSAR processing of the images, including interferogram network selection, coregistration, differential 235 interferogram phase generation and phase unwrapping; ii) deformation time series estimation; and iii) reference point refinement and geocoding.
The main processing steps are described in detail as follows: i) InSAR processing 240 To overcome the decorrelation caused by changes in the permafrost landscape (e.g., freeze-thaw transitions and vegetation phenology), we generated interferograms with only short time intervals. The precision state vectors obtained from the ESA were applied to reduce the effects of inaccurate baselines. All the SLC images were coregistered to the stack reference of 20180807 acquisition for orbit 150 and 20180801 acquisition for orbit 48. After generating a coregistered stack of SLC images, interferograms were generated by each SAR 245 image with its two sequential acquisitions. The temporal baselines of the individual interferograms ranged from 12 to 36 days, and the perpendicular baseline of all the interferometric pairs was < 200 m. Multilooking with 3 pixels in azimuth and 13 pixels in range was then performed to form a square pixel (~40 m) and reduce noise. To remove the topographic phase, The topographic phase was simulated using the SRTM DEM and subtracted from the interferogram. Next, we applied an adaptive spectral filter to produce differential 250 interferograms. To unwrap the differential phase, the SNAPHU Minimum Cost Flow (MCF) phase unwrapping algorithm  was applied. Two orbits covering an area of approximately 450 km×500 km were processed using ISCE software (https://github.com/isce-framework/isce2). To facilitate the processing of data over such a large area, coregistration and conversion between radar coordinates and geometric coordinate systems were accelerated by the aid of a graphics processing unit (GPU) 255 under the Compute Unified Device Architecture (CUDA) framework.

ii) Deformation time series estimation
In this step, the network of unwrapped interferograms was inverted to construct a timeline of line-of-sight (LOS) displacement maps. A multilook 2 by 2 was applied to interferograms to reduce dataset size. Before 260 the network inversion, we applied unwrapping error corrections by bridging reliable regions . The interferograms, generated by each SAR image with its two sequential acquisitions, all have averaged interferometric coherence above 0.7. Appendix Fig. A1 presents an overview of the interferogaminterferogram network and the averaged interferometric coherence. We applied a weighted least square (WLS) estimator to invert the network of interferograms into a time series. During the inversion, 265 interferograms were weighted by the inverse of the phase variance of the whole interferogram . Different from some studies conducted in permafrost environments that presupposed deformation models to help solve the phase time series , we did not preset any deformation and instead obtained the raw phase time series by minimizing the phase residual. The time series of LOS displacements are relative to the first scene of the datasets and spatially relative to the reference point. 270 The reference point was first set by selecting pixels with extremely high temporal coherence greater than 0.99.
After the raw phase time series were obtained, tropospheric delay correction, phase deramping, and topographic residual correction were applied. The tropospheric delay was estimated in the satellite LOS 275 direction using European Centre for Medium-Range Weather Forecasts (ECMWF) Fifth-generation Reanalysis (ERA-5) data. The processing was conducted in PyAPS software . The phase ramps, which might be caused by residual tropospheric and ionospheric delays, were estimated by a 2-D quadratic model based on reliable pixels and removed from the displacement time series at each acquisition.
The systematic topographic phase residuals caused by DEM errors were estimated based on the 280 proportionality with the perpendicular baseline time series . The processing described above was implemented by MintPy  (https://github.com/insarlab/MintPy).
Two indicators evaluated the quality of unwrapped phases and inverted raw phase time series: the phase closure of interferogram triplets and temporal coherence. The phase unwrapping algorithms add an integer 285 number of 2π phase jumps to recover the unwrapped phase. Interferometric phase noise and discontinuities among different coherent regions may lead to the wrong 2π jumps added to the phase field known as unwrapping error. Unwrapping errors can bias the estimated time series. For an interferogram triplet (Δϕ ij , Δϕ jk and Δϕ ik ), unwrapping errors introduce a nonzero integer component Cint ijk in the closure phase C ijk . Therefore, the number of interferogram triplets with nonzero integer ambiguity can be used to detect 290 unwrapping errors: where Δφ ij , Δφ jk and Δφ ik are the three unwrapped interferometric phases generated from the SAR 295 acquisitions at ti, tj and tk, respectively; wrap is an operator that wraps each input number into [− , ); and T is the number of interferogram triplets. A triplet without unwrapping errors has ≡ 0.
The second index, temporal coherence, represents the consistency of the time series with the network of interferograms : where (for N SAR images and M interferograms) Δ is the unwrapped interferometric phase; A is the × ( − 1) design matrix indicating the acquisition pairs used for interferograms generation (consisting of -1, 0 and 1 for each row with -1 for the reference acquisition, 1 for the secondary acquisition and 0 for all other acquisitions ); ̂ denotes the estimated time series; H is an M×1 all-ones column vector; and j is the imaginary unit. 305 Temporal coherence varies from 0 to 1: pixels with values closer to 1 are considered reliable, whereas pixels with values closer to zero are considered unreliable. A threshold of 0.7 is recommended to be used for a dense network of interferograms. In this study, we used a threshold of 0.85; the pixels with temporal coherence below this threshold were masked from the final result.

iii) Reference point refinement and geocoding
In the natural environment, exposed bedrock in flat terrain is normally selected as the reference point; however, because such exposed bedrock is scarce in the study area, we took great care in the selection of reference points. We selected a reference point near the boundary of Selin Co watershed not affected by permafrost in dry and flat terrain. The residual phases from atmospheric distortions and tectonic movements 315 on the TP are difficult to remove completely in such a large-scale extent. However, the effect could be reduced by setting the reference near the study area. The reference area was homogeneous and had interferometric coherence close to 1. Thus, we adjusted the displacements relative to this reference area,; the location of the reference area is marked in Fig. 1. Finally, we geocoded the deformation time series to the WGS84 coordinate system with a grid spacing of 0.0009°×0.0009° and then reprojected it onto the Albers 320 equal area conic system with a 100 m×100 m grid size. To minimize the effects of extreme values, we also applied a 3-size moving window filtering to the deformation time series.

Extraction of periodic (seasonal) amplitudes and inter-annual deformation rates
The surface deformation in permafrost terrain exhibits the characteristics of both long-term linear deformation and seasonal oscillations (uplift in winter and spring and subsidence in summer and autumn) 325 . Therefore, we applied a sinusoidal (seasonal) function plus a linear (interannual) trend to the displacement time series of each image pixel using Eq. (6).
( ) = a ⋅ + ⋅ sin ( 2 ⋅ + ) + c where t is the time interval with respect to the first SAR image acquisition date, a is the long-term deformation rate, A is the periodic (seasonal) amplitude, T is the period of the seasonal undulation (assumed 330 to be one year), is the initial phase, and c is the residual term. a, , and c are the coefficients to be determined. The peak-to-peak periodic (seasonal) deformation is twice the periodic amplitude of . The seasonal amplitude shown in this work refers to the amount of , not the peak-to-peak deformation amount of 2A.

335
For each orbit, we extracted the LOS periodic (seasonal) amplitude and long-term deformation rate pixel by pixel from the deformation time series and then mosaiced the results from the two orbits together. The spatial grids of the incidence angles from the two orbits were mosaiced as well.

Deformation from the LOS direction to the vertical direction
For flat terrain, deformation is caused mainly by freeze-thaw cycles within the permafrost layer and occurs 340 predominantly in the vertical direction. Hence, assuming no horizontal displacement, the observed deformation in the line of the sight (LOS) direction was converted into the vertical direction using Eq. (7) by dividing the deformation value by the cosine of the incidence angle (the incidence angle ranges of orbit 48 and orbit 150 of the S1 satellites are listed in Appendix Table 1A1): 345 where is the observed deformation along the LOS direction, is the deformation in the vertical direction, and is the incidence angle of the sensor.

Conversion from surface deformation into ground ice meltwater contributions
A considerable amount of ground ice is always buried in permafrost regions, especially just below the permafrost table Mackay, 1983;. Thawing of the uppermost 350 permafrost layer is always accompanied by the compaction of sediment and subsidence of the ground surface due to the melting of super-saturated ground ice (French, 2017). Hence, the higher the ice content in permafrost, the larger the surface subsidence occuroccurred as it was thawed. In this study, we assume that the cumulated long-term settlement is equal to the thickness of ground ice melted, and then releasereleased to the hydrological cycle. We converted the deformation rate into the ground ice meltwater via the following 355 steps: i) Masking areas with slope angles >10° On steep slopes, deformation occurs mainly as downslope movements are driven by gravity (Buckel et al., 2021;; thus, we considered only areas with slope angles below 10°. This threshold of 10° was adopted based on previous works Buckel et al., 2021). 360 ii) Masking areas with the LOS periodic (seasonal) amplitudes ≤1.5 mm or long-term velocities ≤2.5 mm/a These thresholds were similarly set based on previous works Buckel et al., 2021) and field surveys. The periodic seasonal threshold was used to mask areas that were unaffected by permafrost activity (because they do not experience periodic frost heave and thaw subsidence) but exhibited continuous surface uplift or subsidence. Areas with periodic (seasonal) amplitudes below 1.5 mm and long-term (inter-365 annual) velocities less than 2.5 mm/a were likely unmoving considering distortions in the InSAR phase and DEM error.
iii) Conversion from surface deformation into ground ice meltwater Selin Co is located in the lowest area of the watershed, and the water released from thawing permafrost eventually enters the lake. The long-term deformation represents the elevation changes produced by ground 370 ice change. Using Eq. (8), by multiplying the grid cell size, we can obtain the ice volume change for each grid cell and then sum among all the cells to estimate the total ice volume contribution. Because the density of water is ~1.0 g/cm 3 and that of ice is ~0.91 g/cm 3 , the ice volume amount is then multiplied by a factor of 0.91 to convert to the corresponding water volume: where ∆ is the change in the ground ice volume,∆ℎ is the elevation change produced by the ice change, is the area of a grid cell, and is the long-term deformation rate in the vertical direction.

Lake water storage changes
We analyzed the changes in the lake water storage of three large lakes within the Selin Co watershed: Selin 380 Co, Ngoin Co, and the combination of Wuru Co & Qiagui Co. Table 5  During this period, the surface water increased at a rate of 0.2 m/a for 4(a-c). Table 3 lists the surface areas 385 of the investigated lakes on the Landsat acquisition dates andSelin Co. Table 6 summarizes the averaged water surface elevation from each lake, considering all the elevation measurements within each year. Selin Co expanded at a rate of 10.3 km 2 /a from 2015 to 2020 and its surface water elevation increased at a rate of ~0.2 m/a from 2018 to 2020.

390
The tracks of the elevation measurements in these years are plotted in Fig. 2(a-c). To calculate the changes in the lake water storage of Selin Co, Eq. (1) was applied taking the areas of 2408.1 km 2 in 2018 and 2441.2 km 2 in 2020 and taking the water surface elevation change of ~0.4 m between these two years; then, the change in lake volume from 2018 to 2020 was estimated, and finally, the annual volume change rate was obtained by dividing the results of these two years. The annual rate of change in the lake volume of Selin Co 395 during 2018-2020 was ~485×10 6 m 3 /a.

for each year. Subfigures (d)-(f) show the ICESat-2 derived elevations. The solid lines indicate the average values of all elevation measurements within the lake on a given date. The light-colored areas show the mean ± one standard deviation. Subfigures (g)-(i) further show histograms of the standard deviations of the water surface elevation of the lake from all acquisition dates. The red dashed line indicates the mean value.
405 Table 63 Lake areas and surface water elevations from . Lake elevations are2018 to 2020, presented as the mean ± one standard deviation of all elevation measurements of the lake within the year. The tracks of the ICESat-2 elevation measurements are shown in Fig. 4 (a)

Deformation at the borehole and GPR survey sites
The deformation time series at the seven borehole sites are shown in Fig. 5. The detected deformation properties are highly consistent with the field surveys. Sites SLC03 and SLC04 displayed the largest seasonal amplitudes of 13.8 mm and 16.6 mm, respectively; among all the borehole sites, these sites are the wettest and have the most developed vegetation cover, i.e., alpine swamp meadows with coverage higher than 90%. 415 These seasonal amplitudes are, to a large extent, caused by the water-ice phase changes in the freeze-thaw cycle and are thus extremely sensitive to the water content in the active layer . Permafrost can hold more water in the active layer than seasonally frozen ground; this is why site SLC05 in seasonally frozen ground displayed the smallest seasonal oscillation of 0.5 mm. Excess ice was found in the SLC04 borehole; correspondingly, the largest subsidence velocity of 18.9 mm/a was observed at this site. Likewise, 420 SLC02 and SLC03, in which ground ice was found, showed different levels of subsidence. Site SLC05, which is seasonally frozen and located in a river valley between two rivers, exhibited an uplift of 3.2 mm/a, which might be because of sedimentation.
We also examined the deformation characteristics at the GPR survey sites. Among the twenty GPR segments 425 where permafrost existence was confirmed, the seasonal amplitudes were all larger than 2.0 mm, with a maximum amplitude of 16.2 mm and an average of 8.2 mm.

Figure 5 LOS deformation time series at the seven borehole sites. Positive values represent uplift, and negative
values represent subsidence relative to the first scene of the S1 datasets. Fig. 6 and showsFig. 7 show the spatial distributions of the periodic (seasonal) amplitude and long-term 435 (inter-annual) deformation rate, both of which are in accordance with the permafrost distribution map (Fig.   1). The boundary in Fig. 1 differentiating continuous permafrost from seasonally frozen ground exactly matches the boundaries observed foron the seasonal amplitudes in Fig. 6(a)(b) and the long-term deformation rates in Fig. 7.6(c)(d). As shown in Fig. 6-7, widespread and large subsidence and large seasonal amplitudes are located in the upstream portion of the Zhajiazangbu subbasin southeast of Mt. Geladandong, where 440 widespread continuous permafrost is present. Appendix Fig. A3A2 provides the deformation map covering the field investigation region.

Figure 7 Map of the inter-annual deformation velocity in the satellite LOS direction. Subfigure (b) shows an
450 enlarged view of the coring area (black sector in subfigure (a)). Dark grey colored "na" means the information could not be retrieved because of no data or decorrelation. Fig. 87 shows a density plot of the seasonal amplitude versus the deformation rate. Subsidence in the Selin Co watershed was normally between 5 and 20 mm/a (see the statistical details in Table 84) but reached 50 455 mm/a in certain regions, reflecting highly excessive volumes of ice and rapid ice loss in this region. The seasonal amplitude ranged between 0 mm and 60 mm within the watershed area. In the Zhajiazangbu subbasin with extensive permafrost, among the areas with deformation rates greater than 2.5 mm/a, 0.1% of them had seasonal amplitudes greater than 30 mm, 2.2% had amplitudes between 20 mm and 30 mm, 24.1% had amplitudes between 10 mm and 20 mm, 23% had amplitudes between 5 mm and 10 mm, and 50.6% had 460 amplitudes of less than 5 mm; overall, the average seasonal amplitude was 6.9 mm.   Fig. 98 shows a map of the potential volume of water contributed by changes in permafrost ground ice. Table   95 lists the statistics of the potential water contribution volume for the Selin Co basin. The main contribution 470 comes from the Zhajiazangbu subbasin, where the permafrost coverage is 66.2%, whereas the contribution from the Boquzangbu subbasin is limited. The potential rate of ground ice meltwater generation was 55.6×10 6 m 3 /a in the Zhajiazangbu subbasin and 1.9×10 6 m 3 /a in the Boquzangbu subbasin, and the overall rate was 57.4×10 6 m 3 /a. Compared to the rate of change in the lake water storage of Selin Co (485×10 6 m 3 /a), we obtained a ground ice meltwater contribution ratio of 11.8%. Converting the potential water contribution 475 volume into the runoff depth within the watershed yields a 2.7 mm/a increase in the runoff depth.  Ratio of lake volume increasement 11.8% Fig. 4(d) clearly shows that the water surface elevation of Selin Co has increased gradually, although seasonal fluctuations are evident as well. Selin Co has 49 ICESat-2 acquisitions during the investigation period ( Fig.   4(g)), and 69.4% of these acquisitions correspond to elevation measurements with a standard deviation smaller than 0.05 m within the lake. This means that the average values of all elevation measurements within 490 the lake could represent the lake's elevation on a given date. To evaluate the accuracy of the resulting lake area extent, we validated the classification results by randomly selecting 1000 pixels from the Landsat images and visually examining the classification results. The classification accuracy reached approximately 0.98.

Uncertainties and accuracies of lake volume changes
Our analysis shows that during the period of 2018-2020, the water level increased at a rate of ~0.2 m/a, the 495 lake area increased at a rate of 10.3 km 2 /a, and the lake water storage increased at a rate of 485×10 6 m 3 /a. These values are compared with those recorded in previous studies in Fig. 10. Lake area change information are from ; water level change information are from Doin et al., 2015;Zhu et al., 2019a); and water 500 volume information are from Treichler et al., 2019;Li et al., 2019). As illustrated in Fig. 10, the expansion of Selin Co was slow before 2000, with lake area and volume increases of 9.1 km 2 /a and 368×10 6 m 3 /a, respectively. Then, in the period of 2000-2005, the lake expanded extremely fast, with the water level increasing at an approximate rate of 1.0 m/a and the lake area and volume increasing at rates of 60.7 km 2 /a and 1576×10 6 m 3 /a, respectively. After 505 2005, however, these rates of increase slowed down, with those of the lake area and lake water storage slowing to only 12.6 km 2 /a and 553×10 6 m 3 /a during 2005-2017, respectively, and the rate of increase in the water surface elevation slowing to 0.2 m/a during 2007-2011 (Doin et al., 2015). Overall, the values retrieved in this study are all within the ranges of the historical values and are closest to the values after the 2010s.

Uncertainties and accuracies of deformation
Two indicators evaluated the quality of unwrapped phases and inverted raw phase time series: the phase 515 closure of interferogram triplets and temporal coherence. Fig. 119 shows the spatial distribution of the number of interferogram triplets with nonzero integer ambiguity (Eq (4)), with the histogram illustrating the distribution of values within the Selin Co watershed excluding glaciers and water bodies. The areas having smaller than three take part 95% of the watershed, while 72.3% of the watershed has value of zero (no wrapping error on all the interferograms). The evaluated the quality of original 520 interferometric unwrapped phases, the unwrapping errors could be further reduced by bridging reliable regions before network revision .

525
bodies. Fig. 1210 shows the spatial distribution of temporal coherence (Eq. (5)), which is used to evaluate the quality of raw phase time series. 99.0% of the watershed has temporal coherence higher than 0.8, 98.1% has temporal coherence higher than 0.85, 96.0% has temporal coherence higher than 0.9 and 89.1% has temporal coherence 530 higher than 0.95.

535
The Selin Co basin was covered by two orbits of S1 images: orbit 48 and orbit 150. Appendix Fig. A2A3 presents the amplitude and long-term velocity retrieved from orbit 48 and orbit 150, respectively. Fig. A4 shows the velocity distributions in the overlapping area. In the large overlapping area covering 270250 km × 5550 km, the results from orbit 150 and orbit 48 reveal identicalalmost the same deformation characters. and value ranges (for long-term velocity -2.6 ± 4.7 mm/a for orbit 48 and -3.3 ± 5.0 mm/a for orbit 150). This 540 confirms that the seasonal amplitudes and deformation rates extracted from the two orbits are robust even though the acquisitions from the two orbits are not the same (identical (Appendix Table 3A1).
In addition, the deformation characteristics were compared with the observations from the boreholes, and excellent agreement was achieved. We also validated the SBAS-InSAR-derived deformation (which has the 545 same processing flow as applied in this study) with the in situ leveling measurements at the Wudaoliang site on the TP; the relative error was found to be 14.8% regarding the long-term deformation rate . Furthermore, we are currently operating an automatic deformation monitoring device in the watershed, and it is expected to provide an independent validation reference.

Uncertainties in the slope angle threshold in the estimation of ground ice meltwater
During the estimation process, we masked regions with slopes steeper than the threshold of 10°. We also tested the impact of setting slope thresholds of 15° and 20° on the results. When using the 15° or 20° threshold, the volume of meltwater released changed by less than 1%. Hence, the setting of the threshold does not greatly influence the final prediction. 555

Uplift displacement signal
In addition to widespread subsidence detected in the upstream of Zhajiazangbu of the continuous permafrost environment, the uplift signal was also observed in the sporadic permafrost environment in the middle stream of the Zhajiazangbu subbasin, also near some drained ponds. Fig. 1311 shows a sector (location marked in 560 red rectangle in Fig. 1,6-7), in which both subsidence and uplift signals are detected. The mean annual air temperature is -2.0°C, calculated based on ERA5-Land air temperature hourly reanalysis data during 2017-2020. For a better interpretation, Landsat optical image and elevation withoverlaid on hillshade were also presented. The spatial distribution pattern of deformation aligns with landscape and topography very well.
Large seasonal amplitude only appears in the vegetated and wet area, which indicates the water storage in 565 the active layer in a certain way. Uplift signals are generally at the slope feet. Fig. 1311(e)-(h) displays the deformation time series of four sites s1-s4 from high elevation to low elevation. Site s2 is stable, viewing from seasonal amplitude and long-term deformation velocity. The subsidence of site s3 is the result of ground ice melting, confirmed by the large periodic seasonal amplitude caused by frost heave and thaw subsidence in the active layer. The uplift signals of site s1 and site 4 are worth exploring, especially site s4. Site s4 has 570 a possibility of being related to ground ice aggradation since it also exhibits moderate seasonal amplitude, but it is more likely related to sediment accumulation or groundwater table rise regarding its location.

580
Previous research normally focuses on the thaw subsidence signal on the TP, and less attention has been given to the uplift signal. Some of the uplift signals might be related to deposition or the rise of the groundwater table. Some of the uplift signals might be caused by ground ice aggradation. A sufficient water supply accompanied by strong evaporation (cooling effect, energy is taken away) might facilitate the upward 585 freezing of previously unfrozen (or seasonally frozen) sediment. Ground ice aggradation is slightly surprising in the overall warming climate of the study area. However, the upward freezing of previously unfrozen (or seasonally frozen) sediment is still possible and may occur because of sediment accretion (e.g., deltaic and alluvial sedimentation) (French, 2017). A previous study  also detected a complex deformation signal in the permafrost on the northeastern TP and hypothesized that the uplift deformation in 590 lowland regions was caused by excess meltwater pooling, which triggered an increase in the segregation of ice near the permafrost table. On the TP, new permafrost forming is detected on the exposed bottom of Zonag Lake . Thus, on the TP, it might be common for degradation and aggradation of ground ice to both occur in permafrost environments, with degradation representing the dominant pattern and aggradation existing in local areas. Currently, most studies focus on permafrost subsidence signals, and few 595 studies have studied permafrost ground ice aggradation and the causes of uplift signals in local environments.
Nevertheless, the uplift signals in the permafrost environment on the TP are worthy of additional research, and further details on the Selin Co basin are expected to be unveiled and supplemented by the next field survey.

Comparison of Selin Co's expansion during 2018-2020 with historical values
Our analysis shows that the lake area expanded at a rate of 10.3 km 2 /a during 2015-2020, the water level increased at a rate of ~0.2 m/a, and the lake water storage increased at a rate of 485×10 6 m 3 /a between 2018 and 2020. These values are compared with those recorded in previous studies in Fig. 12. Lake area change information is from Zhang et al., 605 2020;; water level change information is from Doin et al., 2015;Zhu et al., 2019a); and water volume information is from Treichler et al., 2019;Li et al., 2019). As illustrated in Fig. 12, the expansion of Selin Co was slow before 2000, with lake area and volume increases of 9.1 km 2 /a and 368×10 6 m 3 /a, respectively. Then, in the period 2000-610 2005, the lake expanded extremely fast, with the water level increasing at an approximate rate of 1.0 m/a and the lake area and volume increasing at rates of 60.7 km 2 /a and 1576×10 6 m 3 /a, respectively. After 2005, however, these rates of increase slowed down, with those of the lake area and lake water storage slowing to only 12.6 km 2 /a and 553×10 6 m 3 /a during 2005-2017, respectively, and the rate of increase in the water surface elevation slowing to 0.2 m/a during 2007-2011 (Doin et al., 2015). Overall, the values retrieved in 615 this study are all within the ranges of the historical values and are closest to the values after the 2010s.

Water contribution by permafrost ground ice melting in the Selin Co watershed
In this study, We assumed that the amount of surface settlement corresponds to the release of the same volume of ground ice as a result of compressing the thawing ice-rich permafrost layer and that the released water eventually enters Selin Co; the potential water contribution volume was estimated accordingly. However, the effect of permafrost degradation on the terrestrial water cycle is complicated . A gradually 625 thickening active layer may retain more water in the soil layer; thus, endorheic basins may collect more water during the thawing season. Additionally, the increase in the active layer thickness due to the warming climate could lead to enhanced evaporation. Furthermore, in permafrost terrain, the interactions between groundwater and surface water are restricted. With continued permafrost degradation, the rates of groundwater recharge and discharge are expected to increase since the impermeable barrier provided by permafrost weakens. 630 Ultimately, the melting of ground ice can lead to the infiltration of more surface water into groundwater aquifers in the basin, thereby increasing groundwater storage in the basin .
Based on the detected deformation, we estimated that the volumetric rate of water being released due to the 635 melting of ground ice could reach 57.4×10 6 m 3 /a for the Zhajiazangbu and Boquzangbu subbasins. If the Alizangbu and Zhagenzangbu subbasins arewere included, the volume contributed by ground ice melting could be even larger. According to the permafrost distribution statistics (Table 1) Table 53). Thus, the water supply from other lakes to Selin Co might be limited considering the circumstances stated above and thus was not incorporated into the calculation in the current study. In addition, there are some wet areas facing strong decorrelation (1.9% of the watershed area has temporal decorrelation smaller than 0.85 and masked), where deformation information could not be retrieved. However, these areas usually 650 experience very large subsidence. If taking into account of these areas, the ground ice meltwater contribution could be higher. If only the retrieved subsidence signal is considered, the ground ice meltwater volume contribution is 57.4×10 6 m 3 /a, accounting for 11.8% of the lake volume increase. This result is consistent with a previous model simulation of the water supply of permafrost degradation to the lake volume increase of the endorheic basins on the TP , which estimated the water supply by multiplying the 655 active layer thickening rate and the average ground ice content and revealed that permafrost degradation contributed ~12% of the water supply to the lake volume increase.
Our three-year study period is relatively short compared to the continuous period of the growth of Selin Co from the 1970s to the present. The deformation rate and the lake water storage increasing rate in the other 660 periods might be different from those in 2017-2020. C-band ERS, ENVISAT, and L-band ALOS PALSAR provide historical SAR data from the 1990s to the 2010s and will be processed and analyzed in a future study.

Conclusions
This is the first study to quantify the contribution of ground ice melting to the expansion of Selin Co. We monitored the surface deformation using the SBAS-InSAR technique over a large area of approximately 450 665 km×500 km across the Selin Co basin and then utilized the long-term deformation rates to estimate the potential volume of water being released via ground ice melting. Then, this amount was compared with the lake volume change during the same period, and the contribution ratio was derived. SBAS-InSAR monitoring during 2017-2020 exposed widespread and large subsidence in the upstream section of the Zhajiazangbu subbasin, which is underlain by widespread, continuous permafrost. Subsidence normally occurred between 670 5 and 20 mm/a (78.6%), indicating highly excess ice and rapid ice loss in the region. During the same period, the lake water storage increased at a rate of approximately 485×10 6 m 3 /a, and the potential rate at which water was being released due to ground ice melting reached 57.4×10 6 m 3 /a, contributing 11.8% of the lake volume increase. The uplift signals in some lowland areas with a sufficient water supply and in drained ponds are worth noting, and more details will be revealed through future field surveys. This study is especially 675 helpful forcontributes to explaining the rapid expansion of Selin Co and equilibratingclosing the water balance at the watershed scale. More importantly, the proposed method can be easily extended to other watersheds underlain by permafrost to help us better understand the hydrological changes in these watersheds.

Competing interests 695
No conflict of interest.  "p" denotes the path, "r" denotes the row of the frame, and "D" denotes the acquisition date.