Long-term InSAR and streamflow recession analysis reveal accelerated permafrost degradation in the mining area of the Qilian Mountains

Permafrost underlies about 40 % of the Qinghai-Tibet Plateau (QTP), where climate warming and human activities increasingly threaten fragile alpine ecosystems, necessitating long-term permafrost monitoring. Interferometric Synthetic Aperture Radar (InSAR) enables precise detection of thaw-induced surface deformation, while streamflow recession helps reveal subsurface hydrological changes with permafrost degradation. This study performed a first-time joint analysis of decades-long InSAR surface deformation and streamflow recession to assess the trajectory of permafrost degradation in the source region of the Datong River, an area located in the Qilian Mountains of the northeastern QTP and subject to intensive mining during the 2000s and early 2010s. A data-constrained Small Baseline Subset (SBAS) method was proposed to improve the Sentinel-1 C-band deformation retrievals through integrating a linear–periodic temporal constraint model and using concurrent ALOS-2 retrieved deformation rate as a reference. A consistent long-term (1997–2023) deformation dataset was then generated through combining multi-sensor C- and L-band SAR retrievals. The results reveal minimal surface deformation before the mining, followed by sustained ground subsidence (−15 to $-\mathrm{5}\mathrm{mm}{\mathrm{a}}^{-\mathrm{1}}$) and enhanced seasonal deformation ($\sim \mathrm{20}-\mathrm{60}\mathrm{mm}$) during and after mining, indicating accelerated permafrost degradation. This acceleration coincides with a marked slowdown in the post-mining streamflow recession rate derived from daily discharge data of the upper Datong River, likely driven by thaw-induced increases in basin subsurface water storage and flowpath connectivity. This study provides a first comprehensive assessment of permafrost degradation from both surface and subsurface perspectives, offering valuable insights for integrating remote sensing and hydrological observations to assess permafrost vulnerability.

1 Introduction

Permafrost, a critical component of the cryosphere, plays a vital role in regulating hydrological, ecological, and geomorphological processes in cold regions (Biskaborn et al., 2019; Jin et al., 2021; Zou et al., 2017). The Qinghai-Tibet Plateau (QTP) has the most extensive permafrost coverage in the mid- and low-latitude regions (Zou et al., 2017). Compared with the northern high-latitude permafrost regions, permafrost on the QTP is generally warmer, with ground temperatures exceeding −1 °C in some regions and close to the freezing point, making it highly sensitive to climate change (Mu et al., 2017; Ran et al., 2018). In addition to climate warming, human activities such as mining, overgrazing, and infrastructure development have been exerting growing pressure on the QTP permafrost stability (Hjort et al., 2022; Ran et al., 2022). Under the combined pressure of climate warming and human disturbance, the QTP has undergone widespread permafrost degradation, including active layer deepening, near-surface permafrost loss, and more frequent thermokarst development (Murton, 2021; Yi et al., 2025). To better understand these degradation processes and their environmental impact, long-term monitoring of permafrost dynamics is essential.

Because permafrost is a subsurface feature, it has been challenging to detect its state and evolution using satellite remote sensing (Runge et al., 2022). Nonetheless, permafrost dynamics can be inferred indirectly from surface changes (Zhang et al., 2022). In permafrost regions, the melting of ground ice and the seasonal freeze–thaw cycle of the active layer can lead to long-term surface subsidence as well as periodic ground uplift and subsidence (Liu et al., 2012). Therefore, continuous monitoring of surface deformation can provide critical information for assessing permafrost conditions. Traditional methods, such as leveling and total station surveys, provide highly accurate deformation measurements but are labor-intensive and spatially limited (Tao et al., 2025). The Global Navigation Satellite System (GNSS) offers precise point-based monitoring, but high maintenance costs limit its applicability for large-scale observations (Liu et al., 2025). Recently, Interferometric Synthetic Aperture Radar (InSAR) has become an effective technique for detecting surface deformation in permafrost regions, offering millimeter precision, high spatial resolution, and broad coverage across large and often inaccessible areas (Chen et al., 2022; Liu et al., 2010). However, temporal decorrelation poses a major challenge for InSAR applications in permafrost environments, primarily due to rapid surface changes induced by freeze–thaw cycles, snow cover changes, and vegetation growth (Fan et al., 2025a; Zwieback et al., 2024). This makes it challenging for long-term and continuous monitoring of surface deformation, particularly using shorter wavelengths such as C-band. On the other hand, while long-wavelength SAR data (e.g., L-band) can maintain good coherence over long temporal baselines, the availability of historical satellite L-band data has been sparse, with longer revisit time and less frequent global coverage compared to C-band systems (Wang et al., 2017). Due to the above limitation, the majority of existing InSAR-based permafrost studies have been constrained to relatively short observation periods (<10 years). How to effectively combine the multi-frequency SAR data (e.g., L-band and C-band) to obtain reliable, long-term surface deformation and detect decades-long evolution of permafrost changes is now urgently needed but remains underinvestigated.

Beyond surface deformation, permafrost degradation also profoundly impacts subsurface hydrological processes. As a low-permeability layer, permafrost restricts downward water movement and regulates surface water infiltration and surface-subsurface water exchange (Walvoord and Kurylyk, 2016). With permafrost degradation, enhanced soil water storage and subsurface flowpath connectivity can reshape the subsurface flow production and alter the streamflow recession characteristics at the basin outlet (Bense et al., 2012). Therefore, the long-term streamflow recession process contains valuable information on the basin-wide permafrost state and changes. Previous studies have linked permafrost degradation with increasing winter baseflow, altered ratio of surface to subsurface runoff production, and slowdown in the streamflow recession (Fan et al., 2024; Feng et al., 2022; Walvoord and Kurylyk, 2016; Wang et al., 2018). These studies have revealed profound influences of progressive thawing of ground ice and deepening of the active layer on the river discharge and groundwater changes (Jorgenson et al., 2025). However, few studies have linked the hydrological indicators with InSAR-derived surface deformation to provide an integrated perspective on the permafrost change detection at the regional scale.

The source region of the Datong River, located in the Qilian Mountains of northeastern QTP, had undergone intensive mining activities in the early 2000s and 2010s (Fig. 1; Wang et al., 2022). This area has also experienced significant warming ($\sim \mathrm{0.04}\mathrm{°}\mathrm{C}{\mathrm{a}}^{-\mathrm{1}}$) in the past few decades (1973–2022). To what extent these changes have affected the alpine ecosystem and permafrost state has not been well investigated. Therefore, this study aims to provide a comprehensive assessment of the decades-long (>20 years) permafrost changes in this region from both surface and subsurface perspectives through integrated analysis of long-term InSAR deformation retrievals and the streamflow recession process. Multi-source C- and L-band SAR data were combined to derive a decades-long surface deformation dataset since the late 1990s. We improved the Sentinel-1 C-band InSAR time series inversion through integrating a linear–periodic temporal constraint model and incorporating L-band ALOS-2 deformation rates as a velocity reference to mitigate the influences of a disconnected interferogram network on the C-band deformation retrievals. Long-term daily discharge data at the outlet of the upper Datong River were used to characterize the streamflow recession dynamics and detect subsurface hydrological changes associated with permafrost degradation. By jointly analyzing long-term InSAR deformation and streamflow recession characteristics, this study provides valuable insights into the underlying processes driving permafrost vulnerability.

Figure 1 Overview of the study area. (a) The location of the source region of the Datong River in the northeast Qinghai-Tibet Plateau (QTP) and the in-situ data distribution, overlain on the permafrost map from Zou et al. (2017). The red rectangle indicates the scope of the Muli mining, with the major open-pit mine sites shown as red polygons. The insert panel shows the location of the Datong River basin in the QTP. (b) Elevation map from Copernicus 30 m global Digital Elevation Model (GLO-30, © http://www.opentopography.org, last access: 10 July 2026). (c) Optical image of the Jiangcang open-pit coal mine from Esri | Powered by Esri, with the location indicated by the black rectangle in panel (a). (d–f) Small ponds with a few to ∼30 m in diameter were common in the eastern part of the Jiangcang open-pit mine, as indicated by the UAV images and field photo acquired in July 2025.

2 Study area and in-situ data

The study area is located in the source region of the Datong River, in the Qilian Mountains, the northeastern margin of the QTP. The Datong River flows through a broad, gently undulating valley extending from the northwest toward the southeast of the Qilian Mountains. The source region of the Datong River covers an area of ∼4576 km², and is characterized by widespread (>90 %) permafrost (Fig. 1a). It belongs to the high-altitude permafrost zone of the Altun–Qilian Mountains system on the QTP. The ice-rich permafrost is mainly distributed across gentle slopes in front of mountain ranges (Wang et al., 2020). Permafrost thickness ranges from approximately 50–90 m (Cao et al., 2016), while the active layer thickness varies between 0.9 and 2.5 m (Wang et al., 2020). This area has a mean elevation of 3947±222 m (Fig. 1b) and experiences a continental alpine climate. The mean annual air temperature ranges from −5.25 to −2.19 °C, and annual mean precipitation ranges from 500–700 mm, based on the Third Pole meteorological forcing data (TPMFD, Yang, 2023). Vegetation is well developed, dominated by alpine meadow and alpine swamp meadow ecosystems (Li et al., 2011).

This area encompasses the resource-rich Muli coal field, as shown in Fig. 1a, with an estimated coal reserve of approximately 4.13 billion tons, being the largest coal deposit in Qinghai Province, China (Xiao et al., 2023). It hosts several open-pit coal mine sites (Fig. 1a and c), being the largest on the QTP. This area has undergone more than 10 years of excessive and illegal mining activities, starting around 2003. The completion and operation of the Haergai–Muli Railway in 2009 significantly enhanced coal transportation efficiency, which further accelerated the mining activities (Guo et al., 2023). In response to the severe environmental degradation, ecological restoration efforts were initiated in 2015 by the government. Afterwards, mining activities in the Muli area were gradually reduced and eventually ceased by ∼2020 (Wang et al., 2022).

The in-situ data used in this study include daily discharge and air temperature observations. Long-term daily discharge records covering the warm season (May–October) from 1973–2022 at the Garide hydrological station, the outlet of the Upper Datong River (Fig. 1a), were obtained from the Chinese Hydrological Yearbooks. These data were used to analyze long-term variations in streamflow recession. Without known dams or hydropower facilities in the upstream region, the discharge data at Garide station are minimally influenced by anthropogenic activities, making them a reliable proxy for climate variability and permafrost-related hydrological responses. Daily air temperature data were obtained from two sources: the Qilian County meteorological station (1973–2022), provided by the China Meteorological Administration (CMA), and the adjacent Yakou permafrost site (2016–2022), provided by the Heihe Integrated Observatory Network (Liu et al., 2018). These records were used to characterize the thermal conditions of the study area and to investigate their relationships with surface deformation and variations in the streamflow recession process.

A field survey was conducted in July 2025 to document surface features and changes in the study area using DJI Unmanned Aerial Vehicle (UAV) sensors and historical high-resolution images (Sect. 3.1.3). High-resolution (∼5 cm) multispectral imageries at the Green, Red, Red-Edge and Near-Infrared wavelengths were obtained using the DJI M3M sensor. The UAV images revealed a dense distribution of small ponds with low normalized difference vegetation index (NDVI) and high normalized difference water index (NDWI) values across the flat eastern area of the Jiangcang open-pit coal mine (Figs. 1d–f and S1 in the Supplement) and were used for comparison with the historical optical imagery to identify landform changes. In addition, GNSS interferometric reflectometry (GNSS-IR) data at the Yakou station from Zhang et al. (2021) were used to validate the InSAR-derived deformation time series.

3 Methods

This study leverages multi-frequency InSAR data and long-term stream discharge data to assess the severity and evolution of permafrost degradation in the Datong River source region. Multi-source SAR datasets were first processed to generate differential interferograms. For the Sentinel-1 C-band data, we derived deformation time series using the New Small Baseline Subset (NSBAS) InSAR framework (Doin et al., 2012) with a linear–periodic temporal constraint and an ALOS-2 L-band deformation rate reference to improve the stability of Sentinel-1 C-band inversion. We then generated a long-term InSAR deformation dataset through combining the multi-sensor C- and L-band retrievals to quantify the trends of both interannual deformation rates and thaw-season deformation. Furthermore, we conducted streamflow recession analysis using long-term daily discharge data to characterize variations in the recession time constant, an indicator of thaw-induced subsurface hydrological shift with permafrost degradation. To separate the influence of climate change from mining activities, a multiple linear regression model was developed to predict the climate-driven changes in the recession time constant. The integration of surface deformation retrievals and subsurface hydrological analyses provides an independent, multi-faceted assessment of permafrost degradation.

3.1 InSAR data processing and analysis

3.1.1 Multi-source InSAR data processing

We used multi-source SAR data to derive decades-long (>20 years) surface deformation throughout different mining activity stages (Fig. 2). European Remote Sensing Satellite-2 (ERS-2) C-band SAR data from 1997–1998 captured the pre-mining phase, while C-band Envisat Advanced Synthetic Aperture Radar (ASAR) from 2003–2005 represented the early stage of mining. The intensified mining period was characterized using Advanced Land Observing Satellite (ALOS-1) Phased Array L-band Synthetic Aperture Radar (PALSAR) data from 2007–2010. Subsequent ecological restoration and late post-mining stages were tracked using Sentinel-1 data from 2016–2023 and ALOS-2 PALSAR-2 data from 2019–2022. The characteristics and temporal coverage of the SAR datasets used are presented in Table 1 and Fig. 2, respectively. Due to data availability for the early SAR data and ALOS-2 data in the study area, we prioritized selecting SAR images acquired during the thawing season and SAR images from different years with close acquisition dates to capture the seasonal and inter-annual changes in the surface deformation. For the later period, Sentinel-1 provides dense observations suitable for deformation time series inversion. Although both ascending and descending tracks are available over the study area, the descending track has more frequent acquisitions and shows better temporal continuity, whereas the ascending observations contain several extended temporal gaps (Fig. S2), which would reduce the stability of the time series inversion (Li et al., 2022a). The descending track was therefore selected to ensure better interferogram network connectivity.

Figure 2 SAR data (a) coverage and (b) time span. The optical basemap in panel (a) is from Esri | Powered by Esri.

Table 1 Characteristics of the multi-source SAR data used in this study.

Differential interferograms were generated using the InSAR processing software GAMMA (Werner et al., 2020). The main procedures are shown in Fig. 3. Single Look Complex (SLC) images were co-registered to ensure sub-pixel alignment. All SLC images from different sensors were multilooked to a ground resolution of ∼30 m, consistent with the resolution of the Copernicus 30 m global Digital Elevation Model (GLO-30 DEM, European Space Agency, 2021), which was then used to estimate and remove the topographic phase component. Adaptive filtering was applied to interferograms to reduce phase noise (Goldstein and Werner, 1998). Coherence maps were then computed to assess the quality of each interferogram, and a threshold of 0.3 was applied to mask decorrelated, low-quality areas. Phase unwrapping was performed using the Minimum Cost Flow (MCF) algorithm (Costantini, 1998). Long-wavelength trends, typically arising from orbital ramps and atmospheric contamination, were removed by fitting a low-order polynomial (Yanagiya and Furuya, 2020). The Generic Atmospheric Correction Online Service for InSAR (GACOS) was employed to further mitigate tropospheric delays (Yu et al., 2018). All interferograms were manually checked to ensure high coherence and minimal unwrapping errors over the target area.

Figure 3 Workflow for (a) differential interferogram generation and (b) deformation time series inversion based on the NSBAS-InSAR approach. (c) A comparison of the two data constraint models for the NSBAS framework to mitigate the influence of disconnected Sentinel-1 interferograms during the winter due to snow cover for a selected site (38.0363° N, 99.6344° E).

Different methods were used to retrieve surface deformation time series for different sensors. For SAR datasets with sparse temporal coverage over the study area, including ERS-2, Envisat, and ALOS-1/2, the unwrapped differential interferometric phase was directly converted to line-of-sight (LOS) surface deformation. Details of the used interferograms are provided in Table S1 in the Supplement. For the Sentinel-1 data with much higher temporal resolution, deformation time series were obtained using the NSBAS-InSAR based on the LiCSBAS package (Morishita et al., 2020). LiCSBAS conducts a loop closure check by calculating phase triplets to identify and remove bad interferograms with unwrapping errors (Biggs et al., 2007). A root mean square (RMS) threshold of 1.5 radians was used to reject unreliable interferometric links. After excluding low-quality interferograms, a refined InSAR network (Fig. S3a) with a high regional mean coherence (>0.8, Fig. S3b) was used to derive LOS deformation time series through least-squares inversion. A common reference point (38.0895° N, 99.7243° E) was used for all InSAR datasets to ensure the comparability of deformation estimates derived from different SAR sensors. The reference point was identified from Google Earth imagery as an exposed-bedrock surface outside the main deformation areas and exhibited high coherence (>0.95) and low temporal phase variability across the observation periods.

However, likely due to severe decorrelation caused by wet snow, the Sentinel-1 C-band interferogram network becomes disconnected each winter (Fig. S3a), causing a rank-deficiency problem in the inversion. The NSBAS approach solves this problem by imposing a temporal constraint on the SBAS inversion (Fan et al., 2025b; Morishita et al., 2020). While LiCSBAS defaults to a linear constraint, this is inadequate for permafrost regions dominated by periodic seasonal subsidence and uplift due to seasonal freeze–thaw cycles. We therefore tested two constraint models: a purely linear model Eq. (1) and a combined linear-periodic model Eq. (2) as shown below,

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }u\left(t\right)=c+vt,\text{(2)}& {\displaystyle }& {\displaystyle }u\left(t\right)=c+vt+A_{\mathrm{1}}\sin \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }}{T}}t\right)+A_{\mathrm{2}}\cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }}{T}}t\right),\end{array}$$

where u(t) is the deformation time series, c is a constant offset, v is the linear deformation rate, T is the deformation cycle, generally one year, A₁ and A₂ are the sine and cosine coefficients. The linear-periodic model was ultimately selected as it effectively recovers winter uplift signals within the seasonal data gaps, as shown in Fig. 3c. In contrast, the linear-only constraint misrepresents these gaps as continuous subsidence, introducing significant bias into the estimated deformation rates.

Nevertheless, substantial temporal gaps in the interferometric network during the Sentinel-1 period can still lead to uncertainties in derived deformation rates. To further stabilize the inversion, we developed a data-constrained NSBAS method that incorporates long-term deformation rate from ALOS-2. Because the L-band ALOS-2 data maintain high coherence over long temporal baselines, are less prone to unwrapping errors (Fan et al., 2025a), and provide reliable estimates of long-term deformation rates, the annual deformation rates derived from ALOS-2 were used as a reference for the Sentinel-1 inversion. As the spatial coverage of ALOS-2 is smaller than that of Sentinel-1, we first derived an initial Sentinel-1 deformation rate without applying any velocity reference. A cumulative distribution function (CDF) matching (Fig. S4) was then performed to match the Sentinel-1 rates with those from ALOS-2 over their overlapping areas. The derived correction was extended to the full Sentinel-1 coverage, producing a spatially consistent deformation-rate field that served as a velocity reference in the constrained inversion to obtain the Sentinel-1 deformation time series. The inversion is expressed as:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\displaystyle }& {\displaystyle }\left(\begin{array}{c}\mathit{d}\\ \mathit{\gamma }\cdot v_{\text{ref}}{t}_{\mathrm{1}}\\ \mathit{\gamma }\cdot v_{\text{ref}}{t}_{\mathrm{2}}\\ \mathrm{⋮}\\ \mathrm{⋮}\\ \mathit{\gamma }\cdot v_{\text{ref}}{t}_{N-\mathrm{1}}\end{array}\right)=\\ {\displaystyle }& {\displaystyle }\left(\begin{array}{c}\left(\begin{array}{cc}\mathit{G}& \mathbf{0}\end{array}\right)\\ \left(\begin{array}{cccccccc}\mathrm{1}& \mathrm{0}& \mathrm{\cdots }& \mathrm{\cdots }& \mathrm{0}& -\mathrm{1}& -\sin \left(\frac{\mathrm{2}\mathit{\pi }}{T}{t}_{\mathrm{1}}\right)& -\cos \left(\frac{\mathrm{2}\mathit{\pi }}{T}{t}_{\mathrm{1}}\right)\\ \mathrm{⋮}& \mathrm{\ddots }& \mathrm{\ddots }& \mathrm{\ddots }& \mathrm{⋮}& -\mathrm{1}& -\sin \left(\frac{\mathrm{2}\mathit{\pi }}{T}{t}_{\mathrm{2}}\right)& -\cos \left(\frac{\mathrm{2}\mathit{\pi }}{T}{t}_{\mathrm{2}}\right)\\ \mathrm{1}& \mathrm{\cdots }& \mathrm{1}& \mathrm{\ddots }& \mathrm{⋮}& -\mathrm{1}& \mathrm{⋮}& \mathrm{⋮}\\ \mathrm{⋮}& \mathrm{\cdots }& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{0}& -\mathrm{1}& \mathrm{⋮}& \mathrm{⋮}\\ \mathrm{1}& \mathrm{\cdots }& \mathrm{1}& \mathrm{\cdots }& \mathrm{1}& -\mathrm{1}& -\sin \left(\frac{\mathrm{2}\mathit{\pi }}{T}{t}_{N-\mathrm{1}}\right)& -\cos \left(\frac{\mathrm{2}\mathit{\pi }}{T}{t}_{N-\mathrm{1}}\right)\end{array}\right)\end{array}\right)\\ {\displaystyle }& {\displaystyle }\times \left(\begin{array}{c}\mathit{m}\\ c\\ A_{\mathrm{1}}\\ A_{\mathrm{2}}\end{array}\right),\end{array}\end{array}$$

where $\mathit{d}=[d_{\mathrm{1}},\mathrm{\cdots },d_M{]}^{\mathrm{T}}$ is a stack of M unwrapped interferograms produced from N images acquired at time ($t_{\mathrm{0}},\mathrm{\cdots },t_{N-\mathrm{1}}$), $\mathit{m}=[m_{\mathrm{1}},\mathrm{\cdots },m_{N-\mathrm{1}}{]}^{\mathrm{T}}$ is the incremental displacement vector, G is an $M\times (N-\mathrm{1})$ matrix of zeros and ones that links interferograms to incremental displacements, and V_ref represents the reference deformation rate derived from the ALOS-2 data. γ is the regularization weight controlling the strength of the temporal constraint, which was set to 0.0001 in this study. With such a small value, solutions within the connected parts of the interferogram network are minimally affected by the temporal constraint, which primarily serves to stabilize connections across temporal gaps.

Finally, to reduce residual atmospheric and noise components in the deformation time series, we applied a high-pass temporal filter and a low-pass spatial filter using Gaussian kernels in the deformation time series, with widths of 3 times the average temporal interval and 100 m, respectively.

3.1.2 Developing a long-term surface deformation dataset

Based on the surface deformation time series generated using different SAR sensors, we developed a long-term surface deformation dataset, including the seasonal deformation and the interannual deformation rate. The separation of interannual trends and seasonal deformation was handled differently for different sensors. For Sentinel-1, the complete time series captures the full annual freeze–thaw cycle. Therefore, we derived the annual deformation rate by applying a linear regression to the cumulative deformation time series, while the thaw-season deformation was determined from the detrended deformation time series and calculated as the deformation changes between the onset of thawing (generally early May) and the end of the thawing season (early October).

For other SAR sensors, long-interval interferograms directly measure interannual deformation by comparing similar dates across different years. As the temporal coverage of each interferogram is not always identical, we normalized the deformation to an annual scale by dividing the deformation by the actual acquisition interval (in days) and scaling to one year. For interferograms spanning only a short period within the thawing season, the retrieved surface deformation may not fully represent the total thaw-season signal and thus needs to be scaled to the entire thawing period. However, the traditional scaling approach based on the thawing index (Chen et al., 2023) performs poorly in capturing the actual deformation pattern, as verified using Sentinel-1 data (Fig. S5). For example, scaling deformation from May–July 2019 to represent the whole thawing season significantly overestimated subsidence relative to the observed total subsidence. To obtain more reliable results, we developed a Sentinel-1-based temporal scaling method that directly relates partial-season deformation to the total thaw-season signal. The complete Sentinel-1 time series enables quantification of the proportion of deformation within different sub-periods relative to the total thaw-season deformation. Although the acquisition years differ between the Sentinel-1 and other SAR sensors, we assume the proportion of deformation within a given sub-period relative to the total thaw-season deformation remains stable. This scaling relationship is expressed as:

$$\begin{array}{}\text{(4)}& D_{\text{scaled}}=D_{t_{\mathrm{1}}{t}_{\mathrm{2}}}{\displaystyle \frac{D_{S_{\mathrm{1}}}^{\text{Seasonal}}}{D_{S_{\mathrm{1}}}^{t_{\mathrm{1}}{t}_{\mathrm{2}}}}},\end{array}$$

where D_scaled is the estimated full thaw-season deformation for the other sensor, $D_{t_{\mathrm{1}}{t}_{\mathrm{2}}}$ is deformation measured by the other SAR sensor between t₁ and t₂, $D_{S_{\mathrm{1}}}^{\text{Seasonal}}$ is the seasonal deformation from Sentinel-1, and $D_{S_{\mathrm{1}}}^{t_{\mathrm{1}}{t}_{\mathrm{2}}}$ is deformation from Sentinel-1 over the same period from t₁ to t₂.

3.1.3 Environmental drivers of surface deformation

To explore the environmental controls on the spatial variability of surface deformation, correlation analyses and Geodetector factor analysis (Wang et al., 2016) were performed between InSAR-derived mean seasonal deformation and selected environmental factors, including slope, ground ice content (Sheng, 2020), thawing index, growing-season average NDVI, and thawing season average shortwave infrared (SWIR) albedo. Geodetector is a spatially stratified heterogeneity analysis method to detect associations between a spatial response variable and potential driving factors, and the Geodetector statistic q was used to evaluate the explanatory power of each factor for the spatial heterogeneity of seasonal deformation. The q statistic ranges from 0 to 1, with higher values indicating stronger explanatory power, and is defined as (Wang et al., 2016):

$$\begin{array}{}\text{(5)}& q\left(Y\mathrm{|}X\right)=\mathrm{1}-{\displaystyle \frac{{\sum }_{h=\mathrm{1}}^L{N}_h{\mathit{\sigma }}_h^{\mathrm{2}}}{N{\mathit{\sigma }}^{\mathrm{2}}}},\end{array}$$

where Y denotes the seasonal deformation, X represents an environmental factor and $h=\mathrm{1},\mathrm{\dots },L$ denotes the strata or classes of the variable Y or influencing factor X. N_h and ${\mathit{\sigma }}_h^{\mathrm{2}}$ are the number of spatial units and the variance of Y within stratum h, respectively, whereas N and σ² are the total number of spatial units and the variance of Y over the entire study area.

The slope was derived from the GLO-30 DEM data. NDVI and SWIR albedo (band 7, 2105–2155 nm) were extracted from the MODIS daily Terra NDVI (derived from MOD09GA.006, Vermote and Wolfe, 2025) and MODIS daily albedo products (MCD43A3.061, Schaaf and Wang, 2025). The SWIR band is highly sensitive to surface water because liquid water strongly absorbs radiation in this wavelength range (Wang et al., 2008) and was therefore used as an indicator of surface moisture conditions. The thawing index, defined as the cumulative degree-days above 0 °C during the thawing season, was calculated using the MODIS daily land surface temperature (LST) dataset (MYD11A1.061, Wan et al., 2025) after temporal linear interpolation to fill missing LST data. While slope and ground ice content were treated as relatively static factors, NDVI, SWIR albedo, and the thawing index were derived as multi-year averages over 2016–2022, consistent with the temporal coverage of the Sentinel-1 derived seasonal deformation. All data, including the surface deformation, were aggregated to 1 km resolution before the analysis, and all statistical analyses were restricted to pixels within the study area boundary.

While correlation and Geodetector analyses quantify environmental controls, surface observations help identify localized geomorphic responses to permafrost thaw. Ice-rich permafrost thaw can induce rapid landform changes, often leading to the formation of thermokarst lakes and ponds in lowland areas (Jones et al., 2011; Nitzbon et al., 2020). Field observations in July 2025 revealed numerous small water bodies in the study area (Fig. 1d–f). To assess related surface changes, we analyzed historical and contemporary high-resolution optical imagery. Historical imagery was obtained from a declassified Hexagon Keyhole-9 scene acquired on 27 July 1984 (∼0.9 m resolution) from the United States Geological Survey (USGS) Earth Resources Observation and Science Center. The film data were geometrically registered, radiometrically normalized, mosaicked, and georeferenced using manually selected ground control points in ArcGIS Pro. Contemporary optical imagery was acquired using a DJI M3M sensor, which was orthorectified to produce a high-resolution (∼5 cm) orthoimage for comparison with the historical data.

3.2 Streamflow recession and hydrological analysis

The streamflow recession period, when the discharge declines without new water input, provides a valuable window to investigate potential subsurface hydrological changes induced by permafrost thaw (Cooper et al., 2023). As permafrost degrades, the thickened active layer and the thawing of permafrost enhance infiltration and hydraulic connectivity between supra- and sub-permafrost groundwater, leading to enhanced baseflow contributions to stream discharge (Walvoord and Kurylyk, 2016). Consequently, streamflow in degraded permafrost regions mostly exhibits slower recession rates, reflecting enhanced groundwater storage and prolonged subsurface drainage path compared with those in stable permafrost areas (Fig. 4a). To quantify potential subsurface changes in the study area, we analyzed the streamflow recession dynamics using the recession time constant (Ks), a parameter to describe the normalized streamflow recession rate (Gao et al., 2022). The Ks can be estimated from in-situ stream discharge time series by applying a power-law function that relates the change in discharge ($\mathrm{d}Q/\mathrm{d}t$) to the mean stream discharge (Q) (Kirchner, 2009):

$$\begin{array}{}\text{(6)}& -{\displaystyle \frac{\mathrm{d}Q}{\mathrm{d}t}}=a{Q}^b,\end{array}$$

where a (units: d⁻¹) is the recession rate coefficient associated with the streamflow recession time constant, for which we assume $\text{Ks}=\mathrm{1}/a$. b (dimensionless) is the recession shape index that reflects the shape and nonlinearity of the recession curve. Taking the natural logarithm of both sides linearizes this relationship as:

$$\begin{array}{}\text{(7)}& \ln \left(-{\displaystyle \frac{\mathrm{d}Q}{\mathrm{d}t}}\right)=\ln \left(a\right)+b\ln \left(Q\right),\end{array}$$

A linear regression can be performed on the log-transformed variables to obtain the intercept ln (a) and Ks, as illustrated in Fig. 4c.

Figure 4 Schematic illustration of the link between permafrost hydrology and streamflow recession analysis. (a) Conceptual model of permafrost hydrology and changes in the streamflow recession process with permafrost degradation. (b) Recession period selection based on the daily discharge time series. (c) Linear fitting used to derive the recession parameters.

Using the above methods, we generated a long-term (1973–2022) Ks dataset using the daily stream discharge data at the outlet station of the upper Datong River region (Garide in Fig. 1a). Recession periods were first identified from the daily discharge change series as intervals in which discharge declined continuously for at least seven consecutive days ($\mathrm{d}Q/\mathrm{d}t<\mathrm{0}$). Due to a lack of discharge data during the winter months (November–April) in the upper Datong River basin (due to minimal flow), the analysis was restricted to the period between August and October. Since recession periods identified within a single year may be affected by limited data, hydrological anomalies, or an insufficient number of recession events, they may not reliably represent the streamflow recession behavior. A temporal window of multiple years was generally used for the data fitting. By comparing the Ks time series using 1-, 3-, and 5-year windows (Fig. S6), we adopted a five-year moving-window approach.

Figure 4b shows an example where the data points marked in red represent the selected recession period discharge values. All the discharge data during the recession period within the five-year sliding window are aggregated and used for the linear regression to estimate the mean Ks during this period (Fig. 4c). Since parameters a and b are jointly estimated in the power-law recession model in Eq. (7), a systematic shift in b can alter the magnitude of a and thus influence the interpretation of Ks as a recession timescale proxy. Therefore, the parameter b was also obtained from each regression window to examine its temporal variability and assess its potential impacts on Ks. We further performed a control analysis by refitting the recession model with b fixed at unity (b=1), corresponding to the linear reservoir assumption, to distinguish genuine changes in recession timescale from parameter-coupling effects.

To isolate the possible influence of mining activity from climate change on permafrost degradation, we established a multiple linear regression (MLR) model based on selected environmental variables to predict the streamflow recession time constant Ks under natural climate variation. The model was fitted exclusively using pre-mining data (1973–2002) to ensure that the derived relationships reflect natural climate-driven hydrological variability without anthropogenic influence. The fitted model was then applied to the post-mining period (2003–2022) to predict the Ks expected under climate influence alone. Comparison of these climate-based predictions with Ks values calculated from in-situ streamflow observations allowed us to assess whether the post-mining changes in recession dynamics could be fully explained by climate variability and to identify any remaining anomaly potentially associated with mining activity.

The Ks model predictors included air temperature, precipitation, evapotranspiration, surface net radiation, snow depth, and soil moisture. In situ daily air temperatures were used to calculate the annual mean temperature. The remaining variables were derived from the hourly ERA5-Land dataset (Muñoz-Sabater et al., 2021) and averaged over the study region. Specifically, precipitation and evapotranspiration were summed over the thawing season every year. Surface net radiation was calculated based on incoming solar and outgoing longwave radiation and averaged over the thawing season. Snow depth was averaged over the cold season (November to April). Soil moisture was represented as the thawing-season-average volumetric water content in the 0–7 cm soil layer. All climatic variables were smoothed using a five-year moving average, consistent with the Ks time series.

4 Results

4.1 Spatial variability and long-term trends of surface deformation

4.1.1 Surface deformation retrievals from multi-source SAR data

No persistent or significant deformation within the study area during the early period (1997–2002) before the mining was identified from the deformation maps using ERS-2 and Envisat data (Fig. 5). The ERS-2 data captured ground deformation over two intervals: a short thaw-season period from 28 April–2 June 1997, and a longer one-year interval from 11 August 1997–27 July 1998 (Fig. 5a and c). The Envisat data reveal deformation over a two-year interval, from 18 October 2003–22 October 2005 (Fig. 5e). Although the early C-band SAR data are relatively sparse, favorable coherence was achieved in the flat areas, even with long temporal baselines. However, coherence was poor in mountainous regions and areas with steep terrain (Fig. 5b, d and f). The ability of the C-band data to maintain good coherence over extended periods – up to one to two years – further indicates that surface deformation during this period was relatively minor and stable.

Figure 5 Surface deformation and coherence maps from early-stage C-band SAR acquisitions: ERS-2 (a–d) and Envisat (e–f). Positive deformation represents uplift, and negative deformation represents subsidence.

ALOS-1 L-band SAR retrievals reveal consistent increases in both the seasonal and interannual changes in the surface deformation during the early and peak mining period from 2007–2010 (Fig. 6). ALOS-1 data captured surface deformation during the thawing seasons (July–October) in both 2007 and 2010, as well as cumulative long-term deformation spanning from 2007–2010. The derived interferograms exhibit high coherence across the study area (Fig. S7). The results reveal extensive surface subsidence occurring during both the thawing season and the multi-year period, predominantly concentrated around the open-pit mine and the relatively flat southern part of the basin, indicating permafrost instability likely exacerbated by mining activities. Notably, the thaw-season subsidence observed in 2010 was greater than that in 2007 (Fig. 6a and b), and a cumulative trend of subsidence becomes evident over the study period (Fig. 6c–e). These findings indicate intensified permafrost thaw and active layer deepening.

Figure 6 Deformation maps derived from ALOS-1 L-band SAR data: thaw-season deformation from July–October in 2007 (a) and 2010 (b), respectively; interannual cumulative deformation over 1 to 3 years from 2007–2010 (c–e).

The deformation maps based on the ALOS-2 data indicate an even larger subsidence trend (Fig. 7). In addition to capturing thaw-season and interannual changes of surface deformation, ALOS-2 data also cover the freezing period, specifically from 14 July 2019–22 March 2020 (Fig. 7b). The ALOS-2 derived interferograms also exhibit high coherence (Fig. S8), demonstrating good performance. Consistent with the spatial patterns observed in the ALOS-1 retrievals, the southern part of the study area exhibits widespread and pronounced thaw-season subsidence (Fig. 7a), along with a persistent long-term subsidence trend (Fig. 7c–e), indicating ongoing permafrost degradation. Moreover, the deforming area shows a pronounced surface uplift during the winter season (Fig. 7b), which is consistent with ground heave caused by the refreezing of the active layer.

Figure 7 Deformation maps derived from ALOS-2 L-band SAR data: seasonal deformation dynamics during the thawing and freezing periods (a, b) and interannual cumulative deformation over 1–3 years from 2019–2022 (c–e).

Using the proposed data-constrained NSBAS method, we generated Sentinel-1 deformation time series covering the period from January 2016 to April 2023 (Fig. S9). From these time series, we derived both the interannual deformation rates and the mean thaw-season deformation (Fig. 8a and b). Time-series deformation maps reveal a gradual increase in subsidence beginning in May and a subsequent uplift starting around October, corresponding to the progressive thawing and refreezing of the active layer (Fig. S9). At the two representative sites, P1, with ice-rich permafrost and widespread thermokarst ponds, shows higher subsidence and seasonal amplitude than the drier P2 site (Fig. 8c–e). The duration of subsidence and uplift at both sites follows air temperature changes, with P1 experiencing a longer thaw-season subsidence period (Fig. 8e). A comparison between Sentinel-1 derived deformation and GNSS-IR observations at Yakou station during the snow-free period of 2018 shows good agreement (Fig. S10), confirming the reliability of the Sentinel-1 time-series retrievals. The above results indicate that incorporating L-band ALOS-2 deformation rates as reference information in the inversion model ensured the temporal stability and consistency of the Sentinel-1 time series with ALOS-2 results.

Figure 8 Deformation derived from Sentinel-1 data: deformation rate from 2016–2022 (a); multi-year average deformation during the thawing season (May–October) from 2016–2022 (b); UAV images acquired by the authors in July of 2025 at selected sites: P1 near the Jiangcang open-pit mine with prominent thermokarst features (c) and P2 at the Yakou permafrost station with drier surface and sparse vegetation (d); and the derived cumulative deformation time series for the two sites from Sentinel-1 and ALOS-2 (e).

Spatially, the Sentinel-1 derived deformation rate (Fig. 8a) and the seasonal deformation (Fig. 8b) exhibit consistent patterns with those observed in the ALOS-1 and ALOS-2-based deformation maps. Significant subsidence occurs around the open-pit coal mine and in the adjacent low-lying terrain, where deformation rates range from −15 to $-\mathrm{5}\mathrm{mm}{\mathrm{a}}^{-\mathrm{1}}$. Thaw-season deformation generally ranges from 20–60 mm. Due to the large spatial extent of the map, detailed deformation features of the open-pit mine in the western part are not clearly visible; a zoomed-in view in Fig. S11 shows noticeable subsidence in the adjacent flat terrain. A direct comparison between the L-band ALOS-2 and C-band Sentinel-1 results over their same time period in 2019 further demonstrates strong spatial consistency (Fig. 9). Both maps capture similar spatial patterns of seasonal subsidence, while ALOS-2 shows slightly larger deformation magnitudes and greater spatial variability (Fig. 9a–c). The two maps correlate well (R=0.84; Fig. 9d), with a mean bias of −3.6 mm, a mean absolute difference of 6.5 mm, and a root mean square error (RMSE) of 9.3 mm. Given that the typical accuracy of InSAR-derived deformation is around 5 mm, the observed discrepancies between the two datasets are considered acceptable.

Figure 9 Comparison of surface deformation derived from ALOS-2 L-band and Sentinel-1 C-band data during a similar period of the thawing season in 2019: (a) ALOS-2-derived deformation from 19 May–14 July 2019; (b) Sentinel-1-derived deformation from 19 May–18 July 2019; (c) histogram distributions of the two deformation maps; (d) density scatter plot of the two deformation maps with the correlation and difference statistics, including correlation coefficient (R), mean bias (MB), mean absolute bias (MAB), and RMSE.

4.1.2 Long-term deformation trends derived from multi-source SAR data

The deformation rate maps during different periods show no evident or widespread subsidence within the basin during 1997–1998 and 2003–2005 (Fig. 10a and b). A distinct subsidence trend ($\sim -\mathrm{10}$ to $-\mathrm{5}\mathrm{mm}{\mathrm{a}}^{-\mathrm{1}}$) emerged in the flat area during 2007–2010 (Fig. 10c), marking the onset of notable ground subsidence. Subsidence intensified in subsequent years, ranging from −15 to $-\mathrm{5}\mathrm{mm}{\mathrm{a}}^{-\mathrm{1}}$ during 2016–2023, as revealed by Sentinel-1 data (Fig. 10d). ALOS-2 results during 2019–2022 show a consistent pattern, further confirming the strong and spatially extensive subsidence during the recent period (Fig. 10e). Thaw-season deformation maps show that flat areas initially experienced relatively small seasonal subsidence ranging from about −20 to −10 mm (Fig. 11a), which intensified after 2007 ($\sim -\mathrm{40}$ to −20 mm; Fig. 11b) and further increased by 2010 ($\sim -\mathrm{60}$ to −40 mm; Fig. 11c). Subsequent thaw-season subsidence patterns derived from Sentinel-1 averaged from 2016–2022 (Fig. 11d) and ALOS-2 (2019; Fig. 11e) are broadly comparable to those observed in 2010, indicating the persistence of strong subsidence in the study area. These results suggest that significant permafrost degradation likely began between 2005 and 2007 and accelerated afterwards.

Figure 10 Comparison of average deformation rates normalized to 1 year derived from multi-sensor SAR data across different periods.

Figure 11 Comparison of thaw-season deformation normalized to the full thawing season derived from multi-sensor SAR data across different periods.

4.1.3 Environmental drivers of seasonal deformation

The covariation between seasonal deformation and environmental factors indicates that areas with higher ground ice content, lower SWIR albedo, denser vegetation, gentler slopes, and larger thawing indices generally experience greater subsidence (Fig. 12). Among all the variables, SWIR albedo shows the strongest correlation with seasonal deformation (R=0.65), followed by ground ice content (R=0.60). Consistent with these correlation results, the Geodetector analysis identifies ground ice content (q=0.44) and SWIR albedo (q=0.42) as the dominant factors explaining the spatial variability of seasonal deformation, followed by slope (q=0.40), whereas NDVI (q=0.26) and thawing index (q=0.13) exhibit comparatively weaker explanatory power. The spatial distribution of surface deformation is closely associated with the distribution of surface moisture (indicated by SWIR albedo) and ground-ice conditions (Fig. 13), with significant subsidence generally observed in ice-rich permafrost areas showing densely distributed thermokarst ponds (Fig. 8). In ice-rich permafrost zones, phase transitions between ice and liquid water occur as both the active layer and underlying permafrost thaw, driving strong surface deformation (Chang et al., 2024; Liu et al., 2010).

Figure 12 Covariation of different environmental factors with seasonal surface deformation derived from Sentinel-1: (a) the absolute correlation coefficients and Geodetector statistic q value of different environmental factors with deformation; (b–f): the variation of seasonal deformation with different bins of environmental factors, including thawing index (b); slope (c); Normalized Difference Vegetation Index (NDVI) (d); ground ice content (e), and SWIR albedo (f).

Figure 13 Comparison of the spatial distribution of environmental factors and thaw season deformation: (a) seasonal deformation maps derived using Sentinel-1 SAR data; (b) thawing index; (c) slope; (d) NDVI; (e) ground ice content; (f) shortwave infrared albedo.

Thermokarst ponds or lakes may form in the low-lying areas with permafrost degradation. On the other hand, areas with more surface water have a lower surface albedo and can absorb more heat and continuously conduct heat to the surrounding permafrost, further accelerating permafrost degradation (Yi et al., 2025). This likely explains why areas with greater thawing-season subsidence generally correspond to areas with lower SWIR albedo, with a pronounced north–south contrast across the study area (Fig. 13). This is partly supported by our analysis using high-resolution optical imagery in the study area. For example, at one site near the Jiangcang open-pit mine (i.e., the P1 site shown in Fig. 8), the UAV optical imagery shows a remarkably dense distribution of small ponds with a size generally less than a few meters in the subsiding area. Comparison of the UAV image with the 1984 Keyhole imagery reveals that the presence of water bodies has expanded to some extent (Fig. S12), reflecting the enhanced development of thermokarst landforms in this area, associated with ongoing permafrost degradation.

4.2 Thaw-induced subsurface hydrological changes inferred from streamflow recession analysis

The recession time constant Ks derived from the in-situ discharge data exhibits an overall increasing trend from 1975–2020 (Fig. 14a). Before the mining started, the Ks remained relatively stable, with a slight increase observed during the 1980s and 1990s. However, during the mining period, the Ks increased rapidly and reached the peak around 2013, approximately doubling compared to pre-mining levels, indicating a substantial slowdown in the streamflow recession process. Following the government-initiated ecological restoration efforts starting around 2015, Ks gradually declined but remained higher than pre-mining values. The parameter b exhibited an overall increasing trend from 1975–2020, with the most pronounced increase occurring during 2009–2013 (Fig. S13), coinciding with the period of intensified mining disturbance. This suggests that the recession process likely became more nonlinear (Hinzman et al., 2022). When b was fixed at unity, the temporal variability of a was reduced, suggesting that part of the apparent trend in Ks was affected by the variability in b. Nevertheless, it still shows a similar trend as the original fitting, with a strong decrease in a (i.e., increase in Ks) during the intensified mining period. This confirms that the recession slowdown during this period reflects a substantive change rather than being entirely a mathematical artifact of shifting b.

Figure 14 Long-term variations in the streamflow recession constant (Ks) with the variations in deformation rate and thaw-season deformation during different periods: (a) Ks derived from in-situ streamflow data and predicted from climatic variables along with key climatic variables: air temperature, precipitation and snow depth. The shaded bands indicate the 95 % confidence intervals; (b) correlation between climatic factors and Ks before mining. Climatic factors including precipitation (P), soil moisture (SM), snow depth (SD), evapotranspiration (ET), net radiation (Rn) and air temperature (T); (c) regional mean deformation rate and thaw-season deformation within the overlapping deformation area observed by all sensors across the upper Datong River basin during different periods.

Analysis between Ks and key climatic variables before the mining period (Fig. 14b) reveals that Ks is positively correlated with air temperature, snow depth, soil moisture, and precipitation, and negatively correlated with net radiation and evapotranspiration. Using these climate variables, a multiple linear regression model was constructed for the pre-mining period, which effectively captured the majority of Ks variability (R²=0.73, p<0.001), indicating that Ks changes during this period were predominantly climate-controlled, with air temperature, precipitation, and snow depth showing strong influences. Model predictions during the mining period failed to capture the sharp increase in Ks and the climate-unexplained residual remained positive from 2009–2016, with its 95 % uncertainty interval entirely above zero, coinciding with intensified mining activity (Fig. S14). These results indicate that the observed Ks anomaly during this period cannot be fully explained by climate variability alone and imply an additional mining-related influence on streamflow recession dynamics.

When compared with InSAR-derived deformation (Fig. 14c), Ks shows a coherent evolution with both annual deformation rate and thaw-season subsidence, with correlation coefficients exceeding 0.8 (p<0.01). Before mining activities, seasonal deformation caused by active layer freeze–thaw processes was present ($\sim -\mathrm{19}\pm \mathrm{13}\mathrm{mm}$), but no long-term subsidence trend was detected ($\mathrm{1}\pm \mathrm{4}\mathrm{mm}{\mathrm{a}}^{-\mathrm{1}}$), indicating stable permafrost conditions. During the early stages of mining (2003–2005), there were limited changes in interannual deformation rates ($-\mathrm{1}\pm \mathrm{1}\mathrm{mm}{\mathrm{a}}^{-\mathrm{1}}$). However, as mining intensified, both thaw-season subsidence and long-term subsidence rates increased significantly. The interannual deformation rate reached approximately $-\mathrm{7}\pm \mathrm{3}\mathrm{mm}{\mathrm{a}}^{-\mathrm{1}}$ between 2007 and 2010, with thaw-season deformation measuring $-\mathrm{34}\pm \mathrm{9}\mathrm{mm}$ in 2007 and $-\mathrm{46}\pm \mathrm{10}\mathrm{mm}$ in 2010, which coincided with a sharp increase in Ks, reflecting accelerated permafrost degradation. Ecological restoration after 2015 slowed subsidence acceleration. During 2016–2023, the interannual deformation rate ranged between −11 and $-\mathrm{6}\mathrm{mm}{\mathrm{a}}^{-\mathrm{1}}$, comparable to that of 2007–2010, with a slight increase after 2019. The thaw-season deformation ($-\mathrm{45}\pm \mathrm{10}$ mm) during this period also remained similar to that in 2010. The turning point of Ks occurred around 2013, but due to the lack of satellite deformation observations between 2010 and 2015, it is unclear whether subsidence followed a similar pattern of increase and subsequent decline. Nevertheless, the comparable magnitudes of deformation and Ks before and after this transition indicate an overall good consistency between surface and subsurface responses to permafrost thaw. The persistent subsidence trend after 2019 parallels the recent climate-driven increases in Ks, both responding to continued warming and increased precipitation, suggesting that mining-induced disturbances have amplified permafrost instability and heightened its sensitivity to climate change.

5 Discussion

5.1 Accelerated permafrost degradation during the mining and post-mining period and potential drivers

A decades-long (1997–2023) surface deformation derived from multi-source SAR data reveals accelerated permafrost degradation in the upper Datong River basin, characterized by a significant increase in annual subsidence rates and thaw-season subsidence following the intensification of mining activities. The subsidence rate increased from a stable, pre-mining baseline to approximately −15 to $-\mathrm{5}\mathrm{mm}{\mathrm{a}}^{-\mathrm{1}}$ during the post-mining period. These values are consistent with those reported by Wang et al. (2022), who observed accelerated permafrost degradation around Muli open-pit coal mines, with an average deformation rate of about $-\mathrm{12}\mathrm{mm}{\mathrm{a}}^{-\mathrm{1}}$ in alpine swamp meadows during 2018–2020. Comparable rates have also been documented along the Qinghai–Tibet Engineering Corridor, where anthropogenic disturbances produce deformation of about −20 to $-\mathrm{8}\mathrm{mm}{\mathrm{a}}^{-\mathrm{1}}$ (Lin et al., 2025). In contrast, natural permafrost areas on the northeastern QTP show relatively smaller deformation, averaging $-\mathrm{1.1}\mathrm{mm}{\mathrm{a}}^{-\mathrm{1}}$ during 2003–2011 and $-\mathrm{2.1}\mathrm{mm}{\mathrm{a}}^{-\mathrm{1}}$ during 2014–2019, primarily due to climate warming-induced gradual permafrost thaw (Daout et al., 2020). Overall, deformation rates across the QTP vary considerably depending on surface and subsurface conditions but are typically between −13 and $-\mathrm{1}\mathrm{mm}{\mathrm{a}}^{-\mathrm{1}}$ (Liu et al., 2025). The thaw-season subsidence in the upper Datong River intensified markedly after 2007 and persisted through 2016–2022, with amplitudes ranging from 20–60 mm. These values are substantially higher than those in relatively arid QTP permafrost regions, where seasonal deformation typically remains below 10–20 mm, such as northwestern Tibet and the Hoh Xil area (Daout et al., 2017; Lu et al., 2023), but comparable to the wetter and ice-rich permafrost regions such as the thermokarst-prone area in the source region of the Yellow River (20–50 mm; Li et al., 2023) and the Heihe River Basin (10–60 mm; Peng et al., 2023).

The streamflow recession time constant (Ks) obtained from the long-term discharge data exhibited a rapid increase during the mining period, implying enhanced subsurface water storage and prolonged subsurface flow pathways (Fig. 14a). This response is closely associated with active layer deepening and permafrost thaw (Bense et al., 2012). Previous studies demonstrated that thickening of the active layer can significantly extend the average annual recession duration in the Tahe River Basin, northeastern China (Feng et al., 2022). In-situ active layer observations across four permafrost basins on the QTP demonstrate that supra-permafrost groundwater exhibits seasonal dynamics consistent with the freeze–thaw cycle of the active layer and contributes about 57 %–66 % of total streamflow during the summer period (Qin et al., 2024). Moreover, Evans et al. (2020) reported that in northern Eurasia, enhanced subsurface flow in continuous permafrost regions is primarily associated with active-layer deepening, whereas in areas with discontinuous permafrost, it results from extensive permafrost loss and vertical talik expansion. These processes can be explained by increased infiltration and improved hydraulic connectivity between supra- and sub-permafrost layers as the active layer thickens and permafrost degrades, allowing more water to be retained within the subsurface (Lin et al., 2020; Rogger et al., 2017). A deeper and more connected subsurface flow path also likely reduces the contribution of supra-permafrost groundwater but enhances the contribution of sub-permafrost groundwater to the river discharge, ultimately delaying streamflow recession (Meng et al., 2019). The increase in b during the mining period indicates that the recession process became more nonlinear, further supporting a change in subsurface flow pathways. Permafrost thaw-induced increases in subsurface hydraulic connectivity and drainage heterogeneity may alter the contribution of supra- and sub-permafrost groundwater to baseflow (Hinzman et al., 2022). The interaction of these flow components with different drainage timescales results in a more nonlinear storage–discharge relationship.

The strong coherence between accelerated subsidence and the rapid slowdown in the streamflow recession (i.e., increase in Ks value) during the mining and post-mining period provides compelling evidence that permafrost degradation intensified after the mining activities and is consistent with mining-amplified degradation (Fig. 14c). Before mining, variations in Ks were mainly driven by climatic factors, particularly air temperature, precipitation, and snow depth. Rising air temperatures directly affect the thermal regime of the active layer, thickening the active layer. In addition, thicker snow cover in winter helps maintain warmer soil temperatures, which may promote a thicker active layer in the subsequent thawing season. A thicker active layer enhances the infiltration of meltwater and precipitation, while increased rainfall contributes additional water percolating into deeper soil layers, thereby increasing subsurface water storage and consequently slowing streamflow recession (Jiang et al., 2024; Walvoord and Kurylyk, 2016). Following ecological restoration, Ks began to decline, yet both Ks and surface subsidence remained above pre-mining levels, indicating partial subsurface recovery but persistent thaw influence. Continued increases in temperature and precipitation in recent years may have sustained the slow recession behavior and high subsidence trend. Similar slowdowns in streamflow recession have been observed in many permafrost river basins with recent warming (Jiang et al., 2024; Liu et al., 2024; St. Jacques and Sauchyn, 2009). Moreover, since permafrost degradation had already accelerated after the intensified mining activities, the resulting instability may have increased the sensitivity of permafrost to ongoing climate change.

In addition, areas with greater surface subsidence in the upper Datong River are associated with high ground ice and surface moisture conditions. Based on our field survey, densely distributed thermokarst ponds were widespread in those areas (Fig. 1). Comparison between the historical high-resolution optical imagery and the UAV imagery obtained during our field survey at one of the sites close to the mining area showed the expansion of thermokarst ponds during the past several decades (Fig. S12). Rapid permafrost degradation induced by climate warming and mining activities very likely contributes to the formation and expansion of thermokarst ponds in the low-lying ice-rich permafrost areas, which can absorb more heat and further accelerate permafrost degradation (Yi et al., 2025). This is supported by a close negative correlation between surface SWIR albedo and seasonal subsidence (Fig. 12). However, a further process-based understanding of the link between surface changes and permafrost degradation is still needed. While any surface disturbances, either anthropogenic operations or natural wildfires, can clearly degrade the underlying permafrost, the post-fire degradations seem to be recovered within three to four decades (Cao and Furuya, 2025). It remains to be explored if the degraded permafrost can recover to the original pre-mining state in the future.

5.2 Uncertainties in the data processing and analysis

The long-term deformation results derived from multi-source SAR datasets effectively capture the spatiotemporal deformation patterns associated with mining activities and demonstrate strong consistency with changes in subsurface hydrological processes inferred from in-situ streamflow data. However, certain limitations and uncertainties remain in this study. While multi-source SAR data allow for long-term trend detection, the inherent differences in radar wavelengths (C-band vs. L-band), acquisition geometries (including incidence angles and orbital directions), temporal coverage and observation periods might introduce uncertainties into the derived long-term deformation trends to some extent.

L-band SAR data maintain high coherence over complex terrain and long temporal baselines due to their stronger penetration, producing more reliable deformation retrievals (Abe et al., 2022). In contrast, shorter wavelengths, such as C-band, are more sensitive to surface changes, which can induce temporal decorrelation and unwrapping errors (Fan et al., 2025a). For the C-band Sentinel-1 data, to compensate for interferometric network gaps caused by decorrelation and to enhance the stability of the time-series inversion, we implemented a data-constrained NSBAS approach that incorporated ALOS-2 L-band deformation rates as reference information. This approach takes advantage of the high and spatially continuous coherence of the L-band observations to constrain the long-term deformation trend. Nevertheless, the ALOS-2 reference represents a long-term mean deformation rate, and the linear–periodic temporal constraint may not fully resolve interannual or transient deformation variability within the Sentinel-1 network gaps. The deformation rate derived from the sparse D-InSAR pairs may also be susceptible to local noise and atmospheric effects (Li et al., 2022b). Moreover, ALOS-2 covers only part of the Sentinel-1 study area, and the reference velocity outside the overlapping region was obtained by extending the CDF matching correction, which may not fully represent local deformation variability. An alternative strategy is to connect annual thaw-season networks using interannual Sentinel-1 interferograms acquired near the end of the thawing season (Guan et al., 2024). This approach avoids cross-frequency inconsistencies and requires no additional data. In our study area, these interferograms achieved moderate spatially averaged coherence (0.4–0.7, Fig. S15). Some interferograms, such as the 2019–2020 interferogram (Fig. S16), retained relatively high coherence and spatially continuous deformation signals, demonstrating the potential feasibility of this strategy. However, the coherence maps show strong heterogeneity within most interferograms, and fragmented low-coherence areas can reduce the spatial continuity of reliable phase retrieval and increase susceptibility to phase unwrapping errors (Oliver-Cabrera et al., 2022). For example, localized unwrapping errors were evident even in a relatively high coherence interferogram, appearing as abrupt deformation jumps relative to the surrounding deformation field, as shown in Fig. S16. In contrast, the ALOS-2 L-band interferograms maintained high and spatially continuous coherence (>0.95, Fig. S8) over long temporal baselines, providing a more stable and reliable reference for the long-term deformation trend.

The sensitivity of InSAR-derived phase change to surface deformation is influenced by radar incidence angles and satellite orbit direction. In this study, ERS-2 and Envisat SAR data have relatively small incidence angles (∼23°), whereas ALOS-1/2 and Sentinel-1 data have larger incidence angles (39–40°). Smaller incidence angles generally provide higher sensitivity to vertical deformation (Hu et al., 2024). However, deformation derived from the small-incidence-angle ERS-2 and Envisat data is still substantially smaller than that obtained from the larger-incidence-angle ALOS-1/2 and Sentinel-1 datasets. This suggests that the pronounced subsidence captured by the ALOS-1/2 and Sentinel-1 during the mining and post-mining period is robust and not an artifact of incidence angle differences, supporting the reliability of the observed long-term deformation trends. Assuming negligible horizontal deformation in the study area, we projected the line-of-sight deformation into the vertical direction using the local incidence angles. The resulting vertical deformation trends remain consistent with the line-of-sight deformation patterns, showing significantly greater subsidence during the post-mining period than during the pre-mining stage (Figs. S17 and S18). Moreover, the influence of orbit direction is most pronounced in areas with steep slopes, where horizontal displacements contribute more to line-of-sight measurements. In our study area, the subsiding regions are relatively flat, with slopes generally less than 2°. Therefore, the impact of orbital geometry on the derived deformation changes is considered negligible.

The ERS-2, Envisat, and ALOS-1/2 results were derived from sparse interferograms and therefore represent interval-averaged deformation rates. These observations remain suitable for identifying broad stage-scale deformation differences but provide less information on short-term or nonlinear variations and the precise identification of abrupt deformation changes. In addition, temporal gaps between SAR observation periods leave deformation during these intervals unobserved. The combined results are regarded as a sequence of period-specific estimates rather than a fully continuous deformation record referenced to a common epoch. The incomplete seasonal coverage also introduces uncertainty into the thaw-season deformation estimates. To scale early SAR data to represent the entire thawing season, we assumed that the proportion of deformation occurring within a given sub-period relative to the total thaw-season deformation remains consistent across different years. We evaluated this assumption using the Sentinel-1 time series from 2016–2022 by calculating the standard deviation (STD) and coefficient of variation (CV, defined as the ratio of STD to the mean) of the scaling factor for each sub-period (Fig. S19 and Table S2). Across the four sub-periods, the spatially averaged STD ranged from 0.16–0.23, while the corresponding CV ranged from 0.08–0.11. These relatively low values support the overall validity and robustness of the scaling assumption, particularly within the main subsiding regions. Nevertheless, year-to-year differences in climatic conditions may still affect the relative contribution of deformation during different stages to some extent.

Moreover, our analysis interprets the observed deformation and recession changes as primarily reflecting intensified permafrost degradation after the mining activities. However, mining operations themselves can also cause additional ground movement through confined groundwater leakage or the redistribution of rock mass stresses (Bazaluk et al., 2023), yet such influences remain difficult to quantify in the study area. Additional uncertainties are likely associated with the streamflow recession analyses. While the changes in the recession behavior primarily indicate changes in the subsurface flow production, the lack of winter discharge data restricts a full-year evaluation of subsurface hydrological processes, which dominate in the winter baseflow production. Future studies will extend to other permafrost regions and incorporate continuous year-round streamflow observations to provide a more complete understanding of hydrological responses to permafrost change.

6 Conclusion

This study integrates multi-source InSAR retrievals with long-term streamflow recession analysis to comprehensively assess the decades-long permafrost changes in the source region of the Datong River, an area in the northeastern Qinghai-Tibet Plateau subject to extensive mining in the 2000s and early 2010s. A refined NSBAS-InSAR inversion framework was developed through incorporating a linear–periodic constraint model and ALOS-2 L-band deformation rates as the reference to improve Sentinel-1 C-band time-series retrievals strongly affected by temporal decorrelation during winter. A long-term (1997–2023) deformation dataset was then developed through combining multi-sensor C- and L-band SAR data, including ERS-2, Envisat, ALOS-1/2, and Sentinel-1, to capture both seasonal and interannual variations in the surface deformation. The long-term InSAR deformation retrievals reveal significant subsidence after intensified mining activities, with rates of −15 to $-\mathrm{5}\mathrm{mm}{\mathrm{a}}^{-\mathrm{1}}$ and enhanced thaw-season deformation (−60 to −20 mm), especially in the ice-rich permafrost areas of the upper Datong River basin, indicating accelerated permafrost degradation during the mining and post-mining period. Streamflow recession analysis also demonstrates a pronounced increase in the recession time constant, suggesting recession slowdown following the onset of mining, likely induced by thaw-driven changes in subsurface water storage and flowpath connectivity. The concurrent increase in surface subsidence and streamflow recession time constant during the mining and post-mining periods provides strong evidence of widespread permafrost degradation associated with anthropogenic disturbances. The accelerated permafrost degradation trends revealed by the integrated analysis of InSAR deformation and hydrological data underscore the long-lasting impacts of human activities on permafrost environments and highlight the potential permafrost ecosystem vulnerability under continued strong warming trends in the alpine region.