Volumetric quantifications and dynamics of areas undergoing retrogressive thaw slumping in the Northern Hemisphere

Domain: earth_science
Published: 2025-07-23
Editorial outcome: minor_revision

Abstract
Retrogressive thaw slumping (RTS) is a mass-wasting process characterized by upslope backwasting and rapid thawing of ice-rich permafrost. High-resolution digital elevation models (DEMs) from ArcticDEM enable the volumetric and soil organic carbon quantification of medium to large disturbance areas undergoing RTS (≥10,000 m2) for the Northern Hemisphere. Using DEM time-series analysis and deep learning, we retrieve a total of 2747 disturbance areas undergoing active RTS with a total volume loss of (317.0 ± 0.3) × 106 m3 between 2012 and 2022. Here we show that climatic drivers of RTS activity exhibit latitudinal and regional variations, specifically, the number of precipitation-driven RTS decreases linearly as latitudes increase, whereas temperature-driven RTS increases sharply. Finally, we estimate that 96% of detected RTS thawed ~1.95 × 10–3 Pg carbon per year, equivalent to ~0.2% of annual gradual thaw emission estimates. Our results highlight the complexity of regional RTS dynamics and the importance of high resolution, long-term monitoring efforts.

Introduction
The Northern Hemisphere, particularly the Arctic, is warming at a rate two to four times faster than the global average (refs. 1–3). This rapid warming is driving cascading environmental 4,5 and societal impacts 6, including thawing permafrost 7–9. Permafrost underlies ~15% of exposed land surface of the Northern Hemisphere 10 yet disproportionally consists of 34% of the global soil carbon storage 11,12. When permafrost thaws, the landscape response is largely dependent on ground ice content, ranging from minor elevation change over ice-poor permafrost areas to larger amounts of thaw-subsidence in ice-rich areas causing landscape modifications known as thermokarst 13–16. In recent years, efforts have been made to estimate thermokarst extent and its potential greenhouse gas emissions 17,18, however regional scale empirically-based assessments indicate a wide range of variability in the nature and distribution of thermokarst processes 15. While there is high confidence that climate warming is increasing permafrost temperatures, there is greater uncertainty in the magnitude and distribution of thaw-driven change and estimates of associated greenhouse gas release 12.

Large uncertainties in understanding the consequences of permafrost thaw remain especially for thermokarst on hill-slopes (i.e., thaw-driven mass wasting) due to a lack of large-scale Northern Hemisphere studies concerning the spatial distribution of mass wasting, volumetric rates of subsidence, and a scarcity of information about biogeochemical cycling and emissions of thawed material 19–21. For example, Olefeldt et al. 17 estimated 4.9% of the pan-Arctic (0.91 × 106 km2) to be susceptible to low-to-moderate levels of hillslope thermokarst activity, while Turetsky et al. 18 initiated their hillslope abrupt thaw model spin-ups at 1.4% of all pan-Arctic lands to be affected by hillslope mass wasting (and 3% annual increases of active thaw afterwards). One of the most dramatic and rapidly evolving examples of hillslope thermokarst are areas undergoing retrogressive thaw slumping (RTS) 22–24. Despite increasing interests in understanding the three-dimensional characterization and dynamics of RTS at local to regional scales 8,25–32, the majority of the literature (including recent deep learning studies) focus on one or two-dimensional detection and measurement of RTS (e.g., 33–36) through which thawed volumetrics (and greenhouse gas emissions) cannot be obtained. Quantification of the spatial distribution and rates of subsidence of medium to large RTS (i.e., ≥10,000 m2) is especially important as these disproportionally contribute the most to regional-scale volumetric budgets due to deeper thaw and non-linear area-volume scaling effects 8,31. These uncertainties and information gaps identify an important need for a Northern Hemisphere quantification of the total and annual volume loss of RTS to calibrate future earth system and emission models and validate predictions of previous models.

Thaw slumping describes “the landsliding process involving the melting of ground ice and [subsidence or] collapse of thawed material” (37, p. 91). When this process occurs as backwasting or upslope enlargement it is specified as retrogressive thaw slumping. The permafrost literature describes retrogressive thaw slumping both as a complex set of thermokarst processes as well as a thermokarst landform, termed retrogressive thaw slumps. This identical terminology between process and landform can create inconsistencies in what to map 38 as well as challenge the calibration and validation of machine learning models for large-scale detections 31,38. Henceforth for clarity, we use RTS to denote the process of retrogressive thaw slumping and use the term slumps to denote landforms described as retrogressive thaw slumps 37.

The most common landforms undergoing RTS are thaw slumps; however, RTS can also shape landforms from dead glacial ice in recent ice-cored moraine complexes 39–41. Ice-rich permafrost is required for RTS to occur and as a thaw-driven mass wasting process, the primary process of RTS is the ablation of an ice-rich headwall typically of crescent or elongated shape (Fig. 1c 42;). RTS can be highly dynamic; initializing quickly in a single thaw season and remaining active for decades 33,43, or ceasing activity before potential reactivation occurs 44. However, in the last two decades, intensification of material transfer in response to temperature and precipitation increase has transformed RTS dynamics in many Arctic regions, driving feedbacks that perpetuate retrogressing growth, increase downwasting, and impact downstream environments, leading to the development of areas undergoing RTS to reach between 50,000 m2 and 400,000 m2, creating landforms termed “megaslumps” 8,22. Volumetric enlargement at the retrogressing headwall remains the primary mechanism of thawing ice-rich permafrost, and the disturbance areas undergoing this process are the focus of this study (Fig. 1). A lack of representative RTS training data is often cited as one of the greatest limiting factors to achieve accurate and transferable deep learning models for pan-Arctic detections 45–47, and studies typically overcome this by merging regional datasets from various sources and vintages (refs. 48,49), which might be sparse and biased 50, as well as remote-sensing based and unverified by permafrost experts 38. There is thus a need for a temporally consistent, reproducible, and spatially comprehensive and representative Northern Hemisphere dataset with standardized geomorphic ontological definitions.

Using the ArcticDEM time series 51, here we present the quantification of volumetric change and soil organic carbon release from large disturbance areas undergoing RTS in the Northern Hemisphere for regions ≥60°N (including Alaska, Greenland, Kamchatka Peninsula) for the last decade (2012–2022) and generate a baseline inventory of RTS (i.e., actively retrogressing ice-rich headwalls) with disturbance area and terrain subsidence (total volume loss) estimates as well as annual volumetric yield estimates. These estimates are further correlated with climate variabilities and used to obtain a first-order estimate of soil organic carbon release.

Results
RTS Distribution and Area Estimates
The spatial distribution of medium- to large disturbance areas undergoing RTS (≥10,000 m2; minimum elevation change of 2 m) was obtained using deep learning detection models attuned to thaw-driven geomorphic RTS processes of annual subsidence (see Methods and Fig. 1). Sequential differencing of DEMs from the ArcticDEM time-series (n = 440,949 strips 52;) highlighted incremental surface elevation changes associated with the retrogression of ice-rich headwalls 53, the timing and spatial pattern of which was used as a key distinctive feature during detection (termed the “rainbow” pattern; Figs. 2c, S3e) that distinguished RTS from elevation changes caused by other landscape processes or time-series noise. After manual verification by topical experts, a total of 2747 disturbance areas undergoing RTS were retained across the Northern Hemisphere for which elevation change magnitudes were summarized (Fig. 2a). These verified sites reflect geomorphic change within a 10-year period (2012–2022), characterized by areas of active subsidence through annual thaw-driven headwall ablation and upslope expansion. The median and maximum RTS disturbance area was 20,680 m2 and 249,540 m2 (Fig. 3a), respectively, with a total disturbance area undergoing RTS in the Northern Hemisphere of 78.6 km2 (7860 ha). Compared to the pan-Arctic ground ice map 54, 77% of RTS occurred in areas with high ground ice content (1.37 ha 1000 km−2; Fig. S10; ice content >20%), 5% in areas with medium ground ice content (0.15 ha 1000 km−2; ice content between 10% and 20%) and 18% in areas with low ground ice content (0.20 ha 1000 km−2; ice content <10%). The presence of RTS in low ground ice content areas provides insights on ice-rich permafrost distribution in that localized areas of high ground ice content are still possible within coarsely mapped low ground ice content regions.

The 2747 confirmed disturbance areas were spatially concentrated across several geographic regions (Fig. 2) consistent with results from regional-scale mapping efforts of slumps, such as Banks Island 33, the northwest coast of Victoria Island, Peel Plateau - Mackenzie Delta regions 31,55, Eureka Sound Lowlands 46,56, Novaya Zemlya 57, and western 23 and central Siberian Arctic 30,58,59. In addition to the detection of RTS of common landscape positions and morphologies (i.e., next to fluvial channels or lakeshore- or coastline-based), the deep learning models also identified a series of RTS in Yakutia, Russia (Fig. S13b) that exhibit a different morphology, termed thermoterraces 41, that are morphologically distinct from the typical crescent shape of slumps and are challenging to identify and delineate even by domain experts 38. Other areas undergoing RTS were detected along current glacier margins in northern Greenland, Svalbard, near Klutlan Glacier, Canada; and in Novaya Zemlya, Russia, where glacier dead-ice backwasting occurs (Fig. S13f) (ref. 60). There are considerable geomorphological differences between dead-ice backwasting in recent ice-cored moraine complexes and RTS, yet the small number of large dead-ice backwasting detections (>10,000 m2) were not manually excluded from the inventory as both mass wasting phenomena occur on ice-cored slopes, are thaw related, and exhibit similar elevation change “rainbow” patterns due to their retrogressive, rotational sliding failure mechanisms 39–41. Furthermore, the strong linkage between RTS detections with modern-day glacier margins is compatible with RTS activity across northwestern Canada where terrain encompassing the former glacial margins of the Laurentide Ice Sheet contains relict glacier ice preserved by permafrost 61.

Volumetric Estimates and Growth Dynamics
We estimated terrain subsidence (total volume loss) as a result of headwall ablation and RTS growth over the last decade (2012–2022), with volume estimates consisting of ground ice loss and sediment release that would subsequently flow downslope (e.g., Fig. 1 8,22). The total volume loss for all 2747 areas undergoing RTS was (317.0 ± 0.3) × 106 m3 between 2012 and 2022. The 25 and 100 largest RTS by volume loss (0.9% or 3.6% of total RTS count) contributed 31.0 and 75.3 × 106 m3 (10% and 24%) of the total volume loss, respectively, confirming the disproportionate impacts of the largest areas undergoing RTS 8,31. Regional differences in RTS-driven subsidence indicate heterogeneity in the morphology of RTS, whereby the majority of the largest RTS are distributed in ice-cored moraine of Banks Island and the Peel Plateau, Canada (Fig. 2; and Table 1).

Across the Northern Hemisphere, the median ground surface elevation change for all medium- to large RTS was −3 meters and the median volume loss was 62,484 m3 with a standard deviation of 167,125 m3 (Table 1, and Fig. 3b, c). Most RTS events (86.7%) considered in this study were <200,000 m3 in size (<400,000 m3 95.5%), and average and maximum headwall heights ranged considerably, between 1 to 10 m and 1 to 30 m, respectively (Fig. 3c, d). Larger RTS (by volume) were concentrated within narrow latitudinal bands of 67°N and 81–82°N, while the highest headwalls were concentrated at latitudes 62°N and 65–68°N (Fig. 3e, f). The largest individual RTS by volume loss during the monitoring period occurred 30 km southwest of Aklavik, Canada (68°N) in the Willow River catchment 8, totaling 3,584,000 ± 54,000 m3, a retrogression area of 249,540 m2 (25 ha), and a maximum surface elevation change of −33 ± 1 m. The RTS inventory also addressed knowledge gaps related to the volume loss from thawing headwalls and how that changes with further slump enlargement 8,31. Our study shows a generally consistent (log10) power-law relationship between RTS area and volume loss across different regions in the Northern Hemisphere (Table 1, see “Methods”), with most regions clustered around a scaling exponent of 1.24–1.33. Larger regional contrasts do exist, e.g., scaling exponent ranging from 1.16 (in Plateaus and Plains of Northern Asia) to 1.36 (in Southern Nunavut), indicating a varying degree of non-linearity in the response of RTS expanding volumetrically with increasing area that are likely governed by regional topographic slope configurations and sub-surface conditions.

Temporal metrics provided insights into the annual dynamics of RTS. The majority (75%) of areas undergoing RTS yielded <10,000 m3 annually (<20,000 m3 89.6%), with median yields of 5013 m3 (Fig. 3g). Summarizing annual volumetric yields across regions (named after Provinces and Divisions as defined by 62), RTS located in the Arctic Lowlands (e.g., Banks Island), Mackenzie Mountains, and Brooks Range (incl. Peel Plateau) exhibited the largest yields of 10.0, 4.1, and 4.1 million m3/year, respectively, due to their large population of areas undergoing RTS (Table 1). At a combined 18.2 million m3/year, these three regions contributed up to 56.8% of the observed annual yields across the Northern Hemisphere. The majority of RTS exhibited a relatively stable momentum in their annual yields (i.e., near 0 m3/year2), meaning that annual activity is neither accelerating or decelerating. However, as a whole Northern Hemisphere population, momentum was slightly positive with a mean of 179 m3/year2 (Fig. 3h) indicating a slight acceleration in annual volume loss. These momentum estimates are affected by pronounced annual peaks and lows of activity, but provide useful insights into the net long-term trends of RTS in the Northern Hemisphere as well as sub-regions. For example, across the Northern Hemisphere the decadal trend shows a slight increasing volume loss rate, with total volume loss per year showing peaks in 2016 and 2020 with lows in 2013 and 2018 (Fig. 4b), which are consistent with temperature peaks and troughs in these years (next Section). Temporal trends also varied across regions (Fig. 4c), some characterized by positive momentum with accelerating annual volumetrics, such as Plateaus and Plains of Northern Asia (PPNA) with an acceleration of (0.27 ± 0.05) million m3/year2, Ural Mountains (UM) (UM Division includes Yugorsky Peninsula and Novaya Zemlya Archipelago) with (0.29 ± 0.07) million m3/year2, and Brooks Range (BR) with (0.32 ± 0.08) million m3/year2, while other regions are decelerating (volume loss decreases through time) such as Arctic Lowlands (AL) with (−0.5 ± 0.3) million m3/year2. In addition to regional accelerations and decelerations, some temporal trends were characterized by mid-decadal peaks (e.g., Mackenzie Mountains and Brooks Range) while other regions remained constant (e.g., Mackenzie Lowlands, South Nunavut, and Mountains and Plateaus of Northeast Asia).

RTS Activity in Response to Climate Variabilities
To evaluate the impact of climate variabilities on RTS activity, the correlation between volume loss per year and the air temperature and precipitation time series was calculated for individual RTS (e.g., Fig. 2d), latitudinal and regional groupings, as well as for all RTS in the Northern Hemisphere. Correlations were calculated for the same year as well as the previous year to account for pre-conditioning effects where an extreme summer one year leads to slope failure in the subsequent year (refs. 63,64). Across the Northern Hemisphere a total of 816 RTS (30%) showed strong positive correlations (r ≥ 0.5) (Fig. 4) with same-year air temperature, consistent with current literature 43,56. A total of 329 RTS (12%) showed strong correlations (r >= 0.5) with same-year precipitation trends. Figure 5 highlights latitude-dependent differences in the linkage to temperature and precipitation, with individual RTS binned based on their latitudinal distribution. 71°N marks the inflection point where RTS dynamics were predominantly driven by rising temperatures. The percentage of temperature-driven RTS (influenced by both the same year and previous year temperature) increased sharply from 40% to 88% (as percentage of RTS within their respective distribution bins, Fig. 5) as latitudes increased from 71°N–82°N, and the precipitation-driven RTS decreased linearly from 60% to 13% as latitudes increased from 60°N to 82°N. As shown in Fig. 5, preconditioning effects from climate in the previous year 64 appear prominent at lower latitudes (<65°N), where the percentage of RTS affected by previous-year precipitation (27%) is greater than by same-year precipitation (14%), while at high latitudes (>75°N), the percentage of RTS influenced by previous-year temperature (8%) is significantly lower than those affected by same-year temperature (42%).

The circumpolar assessment in this study demonstrates varied responses between RTS enlargement and climate drivers at the regional scale, agreeing with and further expanding on 35. In the Plateaus and Plains of Northern Asia (PPNA), Northern Mountains (NM), Ural Mountains (UM), Arctic Lowlands (AL, including Banks Island, Victoria Island), 49%, 44%, 43%, 30% of RTS within those regions exhibited strong correlations (r ≥ 0.5) with temperature trends (from the same year); whereas in regions Mountains and Plateaus of Northeast Asia (MPNA), only 16% RTS showed strong correlation with temperature (from the same year). Region Mackenzie Mountains has the highest percentage of RTS (30%) correlated to precipitation trends (from the same year) agreeing with previous studies 22,63, and in region Mackenzie Lowlands, temperature and precipitation played a similar role, with 18% of RTS correlated to temperature and 18% of RTS correlated to precipitation.

Soil Organic Carbon Release
Here we provide the estimate of soil organic carbon (SOC) release from RTS on a Northern Hemisphere scale using empirical volumetric subsidence rates. An estimated 1.95 × 10–3 Pg carbon per year (where 1 petagram = 1 gigatonnes) was released between 2012 and 2022 for 96% of the detected RTS that were within the Northern Circumpolar Soil Carbon Database area (NCSC) 65. The remaining RTS were located outside the database coverage (e.g., Northern Greenland, coastal margins). The NCSC only estimates the upper 3 m of soil, therefore these values are an underestimation for 979 RTS that had mean elevation changes between >3 m and 19 m, however the highest carbon amounts have been shown to generally occur in the upper 3 m and decrease with depth 58.

It could be argued that a large number of small disturbance areas undergoing RTS (area <10,000 m2, small elevation changes) not considered in this study could lead to underestimations of carbon release. To estimate this effect, we used area-frequency distribution modeling and comparisons with detailed manual RTS inventories to estimate the total RTS population size and volumetrics of small areas undergoing RTS (Figs. S11a, S12), resulting in an additional 2244 to 12,514 areas or 18.0–63.8 × 106 m3 (see Methods). The inclusion of small RTS increases the total Northern Hemisphere volume estimate from 317.0 × 106 m3 to a range of 335 to 381 × 106 m3, and shows that the effect of undetected smaller RTS on carbon release is small, comprising only 5–17% of the estimated total volume.

Methods
The ArcticDEM project 51 provides a large collection of very high spatial resolution (2 m), open-access, and time-dependent DEMs covering the entire northern hemisphere landmass above 60°N, as well as all of Alaska, Greenland and the Kamchatka Peninsula, enabling new opportunities for detecting and quantifying RTS. The ArcticDEM data are derived from sub-meter resolution (~0.5 m) stereo-pair images from the Maxar satellite constellation, including Worldview-1 (launched 2007), 2 (2009), 3 (2014), and GeoEye-1 (2008). ArcticDEM DEMs are created using stereophotogrammetry applied to overlapping pairs of images 79, including pairs from the same or different orbital passes (i.e., in-track or cross-track pairs) and combinations of different sensors 79.

The latest release of ArcticDEM (ArcticDEM strips version 4.1, Release October 2023) includes 440,949 time-dependent DEM strips across the Arctic. ArcticDEM data undergo basic automated quality control and filtering for errors due to clouds, shadows, detector saturation 80, water surfaces, and other sources using algorithms developed by the ArcticDEM project team. To better preserve the original data, version 4.1 includes the error masks as separate files instead of applying the masks as in previous versions 52,81.

The overall workflow for RTS detection is summarized in Fig. S1. It consisted of data preparations (Step 1: DEM co-registration and change detection) followed by an initial classification using simple raster-based thresholds (Step 2). This was followed by a multi-step deep-learning approach with iterative model training and validation using several training datasets and data augmentation approaches (Step 3). Finally, a two-step inspection approach using semi-automated decision rules followed by a manual inspection by subject matter experts reduced the sample size to its presented form (Table S4, Step 4).

Step 1: Retrieving surface elevation change from DEM time series
Co-registration of DEMs using bundle adjustment. The DEMs have an internal (pixel-to-pixel) accuracy of ~30 cm (1 sigma) 82. Nevertheless, biases in the horizontal and vertical directions give an external (absolute) accuracy of several meters, caused by the errors in sensor models 83,84. These biases can be removed through coregistration 83,85 and the bundle adjustment of translational offsets among repeat DEMs 82. Validation experiments over unchanging surfaces show that this coregistration reduces the uncertainty in elevation measurements to the random error in the DEMs of ~0.3 m ~0.5 m 82,83,86,87. DEMs with consistently poor pairwise coregistration residuals are then identified as bad data and are removed. Note that to mitigate long-wavelength errors (e.g., tilt) (e.g., 86) in DEMs, we divide these DEMs into 2 km by 2 km tiles and carry out the coregistration and bundle adjustment within the small tiles.

Step 2: Simple classification based on thresholds
Automated change detection algorithm. After coregistering DEMs, we developed automatic algorithms to detect the magnitude and time of the change in the ground surface elevation time series. First, we detect the time of the change by differencing the elevation between sequential observation times. A sudden change in the time series will appear as an outlier/spike in the differenced time series. The advantage of this method is that the differenced time series should cluster around zero (within the noise and seasonal variability) if no abrupt change occurs. The elevation difference is:

y_j − y_j−1 = e_j−1, j=1,2, ...,n−1, (1)

where y_j is the elevation in meters at time t_j, and e_j is the measurement noise. The largest value for y_j - y_j–1 will be considered as a potential event of elevation change if it further conforms to a significance criterion, e.g., larger than 3 times the root mean square errors, and the time range of the event would be between t_j–1 and t_j. The temporal resolution of the event time estimate is mainly limited by the temporal resolution of satellite imagery acquisitions. After estimating the event time, the elevation change can be calculated using a linear model based on the Heaviside step function, as demonstrated in 88. Its uncertainty is estimated based on the posteriori variance of unit weight and error propagation 88. For reader’s convenience, it is reiterated here:

y = a + b × H(t − t_e) (2)

where y is the surface elevation in meters at time t (in years), t_e is the time of the change, H is the Heaviside step function, a represents the surface elevation before the event and b is the to-be-estimated elevation change magnitude, a+b represents the surface elevation after the event.

We then constructed maps of event time of change, elevation change, and uncertainty by applying the above algorithm for each pixel on a regular 2 m grid within a 2 km by 2 km tile. Examples of these maps, shown in Fig. S2 (a, d), display both real surface changes caused by mass wasting and deposits, as well as noise. The elevation change profile (Fig. S2(c)) shows that the magnitude of elevation change increases with the terrain elevation. The detected time of change (Fig. S2(d, f)) demonstrates how this RTS developed retrogressively over the years, where the time of change near the headwall is most recent (2015), and the left of the area had slumping activities in 2013, 2011, and 2010 (earliest where it is furthest away from the headwall). Another characteristic of the time of change map is that the estimated time is locally homogeneous over the depletion zone, but random in other areas.

Fig. S3 demonstrates the sample products from our change detection processing, which include 6 data types for each area. The 6 data types are the elevation change map, the uncertainty of the elevation change, the time of the change T1 (date of the closest DEM data before the detected change), the time of the change T2 (date of the closest DEM data after the change), the number of overlapping DEMs in the area, and the flag of confidence levels for the estimated elevation change at each pixel.

The confidence level is classified for each pixel using the following conditions:
1) Estimated time is homogeneous locally;
2) Elevation change uncertainty < 2 m;
3) Elevation change < −2 m;
4) The elevation change estimates should be larger than 3 fold the uncertainty, and the total number of overlapping DEMs is at least 4;
5) Total area of the selected cluster > 200 m2.

Step 2: Simple classification based on thresholds
From the elevation change maps generated by the change detection algorithm, we then classify potential RTS with a simple thresholding method. The following criteria are used for classifying potential RTS:
1) Surface elevation changes ≤ −1 meter, and total area > 2000 m2.
2) Within the polygon, there is at least one pixel with a high confidence level.
3) The standard deviation of the time of change (T1) within the selected cluster > 1 year.

For two regions under consideration, based on only criterion 1, the total number of RTS candidates is 936,370 for Victoria, Canada (also identified as Region 9) and 9498 for a 50 km by 50 km area in Peel Plateau, Northwest Territories (Table S1). Criterion 3 is designed to identify slumps with multiple thawing periods that display the distinctive “rainbow” pattern. Adding criteria 2 and 3, the total number of RTS candidates is reduced by 8~50 times.

Step 3a: Mask R-CNN for RTS classification with training set 1
To detect RTS in images and generate segmentation masks, we adopt the algorithm Mask R-CNN (Region-Based Convolutional Neural Network) 89, which has been widely recognized as a promising semantic segmentation algorithm across broad scientific applications (e.g., 90–94). Mask R-CNN extends the previous Faster R-CNN architecture 95 to include pixel-level segmentation on top of bounding box recognition. We use the open-source detection and segmentation software Detectron2 96, which is the latest implementation of the Mask R-CNN algorithm 89. Deep learning models in this study are trained with Nvidia A100-SXM Tensor Core GPUs on Google Colab.

Training data preparation
We collected three bands of data as the input for deep learning. As shown in Fig. S4 (a), based on criteria 1–3 only, it is difficult to distinguish RTS from other surface changes, such as ground surface changes around gullies. Hence, two extra bands of information are added as input images. Band 2 is the DEM total curvature 97, which shows distinct patterns near gullies and RTS, e.g., the curvature is high and outlines the headwalls of RTS (Fig. S4g), whereas the curvature shows skeleton features along gullies (Fig. S4c). Band 3 is the time of change map, which shows unique “rainbow” patterns for RTS (Figs. S4h, Fig. S3e), whereas it shows random patterns near gullies (Fig. S4d). The 3-band images (Fig. S4e and i) are then rescaled to integers between 0 and 255, which are the input images for Mask R-CNN.

Each RTS candidate from Step 2 was prepared in a format compatible with the Detectron2 software. We applied a 100 m buffer around the outline, and then the 3-band images, RTS outlines, and bounding boxes were used as the input for Detectron2. The image sample size ranges from ~10,000 m² to 500,000 m². The RTS polygons were labeled as one class ‘slump’. The labeling information including coordinates of polygons outlining the boundary of RTS, bounding box, class types, and file size is exported as JavaScript Object Notation (JSON) files in COCO (Common Objects in Context) format to be read by the Detectron2 software. By scrutinizing candidates from Step 2 in Victoria, Canada, we manually selected the ones with strong “rainbow” patterns as the training and validation data. Training and validation data were randomly split with a ratio of 9:1. The training data contains a total of 900 images, and validation data contains a total of 104 images (Fig. S5).

Training model
Using the training dataset, we fine-tuned a pre-trained Mask R-CNN model (mask_rcnn_X_101_32x8d_FPN_3x.yaml) available in Detectron2’s model library. Referring to the configuration of training parameters in the literature 92,93, the learning rate (cfg.SOLVER.BASE_LR) is set as 0.00025, and the maximum iteration (cfg.SOLVER.MAX_ITER) is set as 48000, and the threshold for prediction confidence score (cfg.MODEL.ROI_HEADS.SCORE_THRESH_TEST) is set as 0.7.

Prediction and evaluation
We adopt the inventory of RTS in the Peel Plateau, Northwest Territories from Segal et al. (2016b) as the ground truth for evaluation. The inventory provides outlines of RTS derived from Quickbird imagery (0.6 m resolution) and SPOT imagery (10 m resolution). We further only selected RTS that show signals from our ArcticDEM data, i.e., the ones that coincide with our RTS candidates from Step 2. The selected 83 RTS in a 50 km by 50 km area in Peel Plateau are used as ground truth (Fig. S5b).

After training, the model prediction is applied to test images for evaluation. The output of prediction for each image is the 2D mask of classified RTS. We remove all clusters with an area of less than 2000 m2. To evaluate the performance of deep learning classification, we adopt the metric intersection over union (IoU). IoU is defined as the ratio of the intersection of two areas to the union of two areas. The threshold for IoU is set as 0.1. A polygon is classified as a true positive if the IoU is larger than the threshold, and a false positive if otherwise. As shown in Fig. S5b and Table S2, there are 83 true positives in this 50 km by 50 km tile area in the Peel Plateau. The criteria 1–3 in Step 2 can successfully recover 77 true positives yet failed to filter out 1161 false positives. The deep learning results can reduce the number of false positives by 23%. We further combine the two solutions, by only selecting candidates that have more than 50% of areas classified as RTS by deep learning. This combination of deep learning and criteria 1–3 reduces the number of false positives by 3 times from 1161 to 357. This combination, by requiring the agreement of two methods to be better than 50%, removes a large number of false positives. As demonstrated in Figure S6a, for the ground truth Slump FM327, the cluster from criteria 1–3 and deep learning show a good agreement, yielding a true positive. For the uppermost cluster, the deep learning indicated pixels (red) are much less than that from criteria 1–3 (white), which indicates bad agreements, leading to a rejection of that cluster. However, a substantial number of false positives still remain, leading to a low precision of 17% (Table S2).

Step 3b: Mask R-CNN for RTS classification with training set 2
After the first iteration of applying the deep learning algorithm (from training dataset 1) on several regions, we further manually selected a group of RTS across broader regions as a second training data. The manual selection of training data sets was conducted in QGIS, by visually scrutinizing each candidate shapefile on top of the time of change base map. The candidate shapefile that encloses an area with the “rainbow” pattern is classified as a true RTS. The dataset is randomly split into training and validation sets at a 9:1 ratio, resulting in 1660 training images and 189 validation images (Fig. S7).

Data augmentation and ensemble modeling
Considering that the available training data size is relatively small, data augmentation techniques (e.g., 91,92) were applied to minimize overfitting. The adopted data augmentation includes rotation, translation, rescaling, and flipping, which increases the sample size by 25 times. The precision and recall of five different models are compared for Peel Plateau as shown in Table S3. The overall issue is that the results have too many false positives, which yield a low precision. Data augmentation in fact unfavorably increases the number of false positives, e.g., 524 vs 357 for model 1 aug (with augmentation) and model 1 (without augmentation). Using the training data set 2, the number of true positives increases slightly, but the number of false positives also increases. The ensemble model is retrieved by applying a logical AND operator to all individual model outputs. Since individual models may detect different false positives, when combined, the overlap of common false positives reduces, hence enhancing the precision (Fig. S6b). For the test area in Peel Plateau, the ensemble combination of five models (Table S3) reduces the number of false positives from 357 to 210. Yet, the precision remained low, only 24%.

The prediction performance shows regional variations due to different land cover and environmental settings such as vegetation cover, snow and ice cover, and water bodies. Further evaluation is carried out in Eureka, Canada (Fig. S5c 56), which represents the Arctic desert environment with little vegetation. There are 45 RTS as ground truth in this test area manually examined by subject experts. Due to minimal impact from vegetation changes, the ensemble model achieves a considerably higher precision of 65% and a recall of 77% (Table S3). The performance differences across land cover types indicate limited spatial transferability of our method.

Step 4
After deep learning, we carried out a full visual inspection to ensure reliable results and achieve high precision. We visually assessed the polygons based on Google satellite images and the “rainbow” patterns in the time of change maps derived in this study. Additional visual inspection using ESRI satellite imagery base map, Yandex Maps satellite imagery base map and multi-year PlanetScope satellite imagery were carried out. Randolph glacier inventory maps 98 were also used to mask out glacier areas. We further used a priori information for references including ground ice content maps 54 and published RTS (Fig. S8). The final selection of RTS is listed in Table S4, which includes a total of 2747 RTS.

Volume estimation
We divide the entire RTS into sub-areas by the time of change in years. For each year, we estimate the volume change per year using this equation:

volume per year = (Σ A_i * h_i) / Δt

where A is the area of a single pixel (e.g., 4 m2), h_i is the height change at each pixel, and N is the total number of pixels within the RTS of a given year. Δt is the time span in years. It’s noted that satellite sampling interval may undesirably scale the volume loss rate. For example, if there are two acquisitions in April and September, the satellite interval will yield a half-year time difference, where it actually means the change occurred in one summer year. Hence we use the number of summer years as the time difference instead of using the precise satellite acquisition dates. Here, Δt represents the count of summers (summer defined as the period from June 20th to September 22nd each year) within the measurements’ time span. If an acquisition occurs before the summer of a year, the time is labeled as the previous year.

Following the method in 80, the uncertainty of volume per year is calculated as:

σ_v^2 = (Σ A^2 * Σ N_ij * C_ij) / Δt^2

Where σ_v is the standard deviation of the estimated volume per year. C_ij is the covariance between height error at pixel i and pixel j. The covariance can be empirically calculated using the periodogram method 80,99 (Fig. S9).

Population size and volumetric estimation of small RTS
Two avenues were pursued to model the total volumetrics of small RTS (<10,000 m2), namely area-frequency distribution modeling and through comparisons with detailed manual RTS inventories.

The active RTS detected in this study yield a power law area-frequency distribution (Fig. S11a), which can be modeled using equation:

log_10(p) = α − βA

where A is the RTS headwall expansion area observed between 2012–2022 (m2), p is the probability density, α and β are coefficients. The fitting of the area-frequency distribution yields α=2.5 ± 0.2, β=(2.7 ± 0.5) × 10−5. This power law distribution indicates that the number of RTS increases as the headwall expansion area decreases. The modeled parameters predict that the number of small active RTS (headwall expansion area <10,000 m2) is about 2244 (upper-bound: n=3444). Considering the area-volume relationship (next section), the estimated total volume of small active RTS is about 18 × 106 m3 (upper-bound: 27 × 106 m3), which is a small percentage (5–8%) of the estimated total volume of small and medium to large RTS ((317.0 ± 0.3) × 106 m3). This prediction is consistent with the observation of 83 RTS in Peel Plateau (Fig. S5b), where small RTS (2000 m2 ≤ expansion area < 10,000 m2) represent 45% of the total count, yet only 12% of the total scar area.

The area distribution obtained in this study was compared to eight detailed RTS inventories. To do so required careful assessment of RTS ontologies and definitions used in the respective inventories. For example, 30’s inventory represents annual average area change using medium resolution DEM differencing, whereas 31’s multi-temporal inventory leveraging high resolution data was used to manually calculate annual average active area changes (n=186) given that this inventory method also included inactive and partially active RTS. Reference 100’s inventory reflect 2002 and 2018 epochs, the areas of which were differenced to obtain annual average change estimates. Furthermore, other inventories commonly cited in the literature, such as 24,45,55,63,101,102 are manual single-epoch delineations that typically reflect larger RTS scar zone areas, and not solely the area of active subsidence measured within a specific time period as in this work. We multiplied the scar area estimates with a mean standardized growth rate of 1.34% (i.e., twice the long-term rate observed between 1950–1960 and 2004–2008 in 55 to account for RTS intensifications in the last two decades) to obtain an estimate of annual area change. Then all samples of the eight inventories were prorated to reflect a 10-year change estimate consistent with this work. The eight inventories (Fig. S12a, b) were characterized by a substantial number of small RTS (with active subsidence area <10,000 m2 over a decade), specifically 82% of total counts (n=6736 small RTS versus 1520 medium-to-large RTS). This estimate does not include inactive RTS or those active for less than a decade 35. These percentages differed among source imagery resolutions and among studies (Fig. S12c, d). If these percentages are generally consistent across the Northern Hemisphere we would expect that the number of areas undergoing RTS detected in this study (n=2747) represents about 18% of all active RTS over a decadal period, for a total Northern Hemisphere population of small RTS of 12,514. Considering the area-volume relationship, the estimated total volume of small active RTS based on this population estimate is 63.6 × 106 m3, which is 17% of the combined estimated total volume of small and medium to large RTS.

Area-volume relationship
The volume of thawed material scales exponentially with RTS area 8,31. Rates of RTS expansion and intensification across scales (or landslides in general) are typically quantified using area-volume allometry, which is a structural trait and an important indicator of potential nonlinearity effects 31,103. The 2747 active RTS from this study provided a precise estimation and interpretation of allometric model coefficients for regional populations as well as the Arctic as a whole using high spatial resolution data:

log_10(V) = a + δlog_10(A)

where, A is the retrogressing headwall area (m2), V is the corresponding volume loss (m3), a is the scaling factor, δ is the scaling exponent 31. As shown in Fig. S11b, the power-law relationship between headwall area and volume is evident and convergent in our product. A least-squares adjustment yields δ=1.30 ± 0.01 and a = −0.83 ± 0.05, with most regions clustered around δ=1.24 to 1.33 (Table 1) indicating a consistent and non-linear headwall yield with slump enlargement across most of the Arctic. RTS belonging to the Plateaus and Plains of Northern Asia and Southern Nunavut represented allometric endmembers, with δ=1.16 and 1.36, respectively (Table 1). In the former region δ suggests that headwall volume loss grows primarily with increasing planform area with little increase in the absolute depth (i.e., nearing linear growth model 31). In the latter region δ suggests that the relative depth of thaw grows non-linearly but not quite in proportion to the planform area (i.e., δ=1.5; self-similarity growth model 31). Small differences in the allometric scaling exponent represent large volume differences, e.g., applying δ=1.5 instead of δ=1.4 produces a volume difference of a factor of 2 104,105. This sensitivity also means that δ is difficult to directly compare against other studies using a landform-type definition of RTS, with 8 and 31 estimating δ of 1.45 and 1.36, respectively, based on very high resolution data in the Beaufort Delta (Canada), or using coarser data across several regions indicating δ ranges between 1.05 to 1.37 30.

Climate variability data
To estimate Arctic air temperature, we used the surface air temperature (at 2-meter height) estimates from the hourly ECMWF Reanalysis version 5 (ERA5) data with 0.25° grid size 106. CRUTS (Climatic Research Unit gridded Time Series, v. 4.07) surface air temperature dataset 107 was also used for comparison (figure not shown). The precipitation data were from the hourly ERA5 dataset (total precipitation) 106. The total precipitation is the accumulated liquid and frozen water, comprising rain and snow, that falls to the Earth’s surface. The mean values are calculated for all RTS sites during the selected summer months (July and August).

Soil organic carbon estimate
We used the Northern Circumpolar Soil Carbon Database (NCSCD) 65,108 to estimate the amount of soil organic carbon (SOC) released in the upper 3 m of permafrost by RTS. We used the shape file version of the dataset and the NCSCD provides a SOC mass (SOCM) estimate for each polygon within the database. We assume that the proportion of SOCM to be constant within each dataset polygon, therefore we calculated the area proportion for each individual RTS and applied this value to SOCM values. The NCSCD divides the SOCM for each 1 m interval for the upper 3 m of soil. Based on the mean elevation of individual RTS, we summed the SOCM data layers accordingly, meaning, for example, RTS with a mean elevation change of 2 m only considered the SOCM amount in the upper 2 m. Due to data availability, RTS with elevation changes of greater than 3 m still only considered the SOCM amounts in the upper 3 m. A total of 10 RTS had mean elevation changes of 1 m, 857 RTS of 2 m, 793 RTS of 3 m and 979 RTS of greater than 3 m, respectively. There were 108 RTS that were located outside of the NCSCD dataset and were not included in the SOCM estimate.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript by Dai et al. is a highly ambitious piece of work that presents novel information on the dynamics of RTSs at the Pan-Arctic scale. By using large amounts of ArcticDEM data, the analysis covers the entire Pan-Arctic at a very high spatio-temporal resolution. The approach shares aspects (using ArcticDEM at Pan-Arctic scale with deep learning approaches) with the work by Huang et al. (2023 ISPRS Journal of Photogrammetry and Remote Sensing, https://doi.org/10.1016/j.isprsjprs.2023.10.008). However, the present manuscript goes beyond identifying RTSs by presenting spatially and temporally detailed volumetric estimates of soil loss and associated carbon release. The found magnitude of carbon release from RTSs provides an important new viewpoint to the globally relevant issue of permafrost degradation and associated carbon budgets. The manuscript also makes an interesting point by showing how volumetric RTS activity demonstrates latitudinal and regional differences in their responses to climate, but also that ecosystem properties markedly complicate this connection.\n\nA comment that needs to be made is regarding the relatively weak classification precision of the presented deep-learning method. High number of false positives is typical for approaches using automated recognition of landforms, and I appreciate the authors for addressing this challenge and pursuing to improve the automated mapping by combining solutions (change detection algorithm and deep learning) and using additional training data with data augmentation methods. Yet, also after this model precision remains low and coming up with reliable results requires undoubtedly laborious visual inspection. On the other hand, the outcomes of the visual inspection are also affected by which locations are selected as candidate RTSs with the initial model and augmented models. Please see my specific comments regarding this aspect below.\n\nThe results are interesting and novel and should be of broad interest to the permafrost research community. The manuscript has an easy-to-follow structure even though the methodology presented involves several processing steps. The figures are clear and interpretable, and the language seems to me (as a non-native speaker) to be in very good order. Below I raise a couple issues and add minor suggestions for corrections and improvements. I recommend publication of the article after my raised issues have been addressed.\n\nSpecific comments:\nrr166-183: Here, the authors describe how the mapped RTS distributions are consistent with previous studies. I would also like to see a discussion of whether there are any regions or environments where RTS activity could not be mapped with a similar precision, or if the results contradict previous findings. Now, the discussion of uncertainties centers on the large number of false positives but the potentially missed RTS activity is less addressed. Although, the limitations in detecting some established activity at megaslumps due to narrow horizontal scale of activity was satisfactorily explained, which indicates that the approach may not have high recall in all conditions. Do we know if model evaluation statistics (from prediction and evaluation with Victoria-Peel Plateau data) hold for the Pan-Arctic area or can there be regional/environmental variation?\nr517: The split ratio of 9:1 sounds reasonable. Were other ratios tried in model training to examine their effect on overall accuracy?\nrr563-564: Were the samples of 1660 training and 189 validation images also here split randomly based on the manually selected and inspected 1660+189 images? Please elaborate because models should be affected by how well the variation in the modelling data is represented.\nrr576-578: I find it quite surprising how the combination of five models (original + augmented) increase precision so drastically. Model ensembles often do outperform single models, but the improvement is surprisingly large. Maybe single models identify different FPs and when put together, there are not so many common FPs, and thus the precision is higher? Could this be elaborated?\nr745: I could not access the produced data via the (supposedly) DOI link (Mendeley?). Are the data available already?\nFig. S7-8. The final detected RTSs are seemingly concentrated in the training data (both 1 and 2) locations (as should be expected) while relatively little RTSs are detected in some areas known to have abundant occurrences. For example, only 13 RTSs have been detected in Alaska. I was wondering why there were no training images in Step 3b selected from Alaska at all? Could this be because the initial model from Step 3a could not predict potential RTSs there due to different environmental conditions in those areas, and consequently that also Step 3b results in very little RTSs in the western part of Brooks Range, for example? Or are these possibly “missing” RTSs simply filtered out in the classification in Step 2 for being too small or not active during the study period?\n\nTechnical touches:\nr99. Replace eg with e.g.\nFig. 1 a, b. Imagery source not visible, please enhance.\nr223: 3gh -> 3g,h\nTable 1: The pink and grey colors are hardly separable. Consider adding contrast.\nFig. 4: The legend for the correlation coefficients has different ranges to those presented in the figure caption and elsewhere in the text. Which are correct?\nrr355-360: This sentence could be more clearly formulated. Is there a comma missing? Consider revising.\nr643: Please elaborate the idea behind selecting the growth rate of double the long-term rate in Segal et al. 2016?\nr652: “…for a pan-Arctic total of ±12,514”. Please check if this is correct. Are confidence intervals/standard deviation for the value missing here or how should one read the value?\nr940: Please check the text font.\nFigs. S7: This is a bit rough figure. Consider using the same longitude of origin for the polar projection that was used in the other maps for better comparability. Please, also acknowledge the imagery source and consider adding a legend in which training and validation images are preferably presented with unique colors.", "reviewer_2": "The manuscript by Dai et al quantify retrogressive thaw slumps (RTS) in the northern permafrost area in terms of location, time of development and volume. The data is generated by comparing multitemporal DEMs based on ArcticDEM and identifying RTS based on a machine-learning algorithm. The data set not only identified polygons of RTS, but also the volume of the landform. This in turn could be used to estimate possible carbon release due to permafrost thaw and erosion. This combination makes the manuscript and its data base highly valuable for the research community, and this study is as far as we know unique and absolute worth publishing, even if the time series of 10 years is too short to evaluate real climate trends.\n\nWe are not experts on the machine-learning algorithms and let this part be evaluated by other reviewers. However, the following major comments are given and may be considered:\n• Title: Be precise, you deal (at least according to your writing) almost exclusively with retrogressive thaw slumps (RTS), why write then “large permafrost-thaw driven mass wasting events”? See also comment below.\n• Type of RTS: We are not always sure what landforms you have detected. Of course, the big RTS in Arctic lowlands are captured and exemplified in images in the paper. However, in the Arctic islands, e.g. on Svalbard, you mentioned large RTS developed (by thermal erosion?) in large recent (< 250 years old) ice-cored moraine complexes in front of recent glacier fronts. Those those would have a quite different development history than ground-ice induced thaw slumps, even in areas of buried old Pleistocene glacier ice (there is a large literature base on this thermal erosion in ice-cored moraines on e.g. Svalbard). Maybe, you should show different types of identified RTS from the various regions? We wanted to have a look in your data for e.g. Svalbard, but the link provided failed.\n• Correlation of events to temperature and precipitation: First of all, the results of the correlation analysis, which is interesting, is in the Discussion (chapter 3.1.), but it should be a Result. Consider moving the first part of that chapter (l. 252. onward). You then write that 60% of the RTS are neither correlated with trends in ERA temperature or precipitation. I am not sure how you made this correlation, but an event one year is not necessarily correlated with climate the same year, but with the weather development the previous year. E.g. an extreme summer one year may lead to weakening of the structure and failure in one of the subsequent years. You should consider this in your discussion (se pre-conditioning papers in high-mountain permafrost slopes, e.g. here https://iopscience.iop.org/article/10.1088/1748-9326/ad3c55). This pre-conditioning may be more important further south than in the northernmost parts?\n• Soil organics release: These are nice data, again, why do you not present these as a result, and then discuss them later.\nLine 341 ff: The reference Tanski et al. (2017) does not support the statement of less SOC release. Here, the numbers are just small because SOC is buried again. However, in this study, the deposition of SOC is not taken into account, so that your numbers should be compared only to the released SOC.\n• Discussion in general: The Discussion in general could benefit from a more detailed elaboration of the results from Fig. 3, which is nice, but not discussed further in the paper (just presented as a result). E.g., why do you observe such areas, volumes, depths, rates and accelerations? Is air temperature and precipitation the only explanation as you focus on this in the discussion? What about ground ice contents, topography, geomorphology etc.?\n• Comment on the so-called “Rainbow pattern”: What about RTS where activity started or ended during the observation period and thus show no or only a weak rainbow pattern? Are these excluded? If that is the case it should be discussed in the uncertainties.\n• Uncertainties: A major result of the discussion is the relationship between RTS activity and climate, and the influence of RTS on carbon release. This discussion is based on ERA5 data and SOC map by Hugelius et al. Here, we miss a more in-depth description of the uncertainties in ERA5 and the SOC map, which of course influence the conclusions. soil organics carbon map.\n• Figures: Some of the figures are a bit messy and should be revised. A list of suggestions are given below in the minor comments.\nMinor comments are as follows:\nFigures and tables:\nFig. 1: Scales are lacking in a, c and d. Lat/long are missing in at least a and b.\nIs c and d the same RTS? If yes, which one of the RTS shown in b?\nFig. 2: Scale is missing in a, no lat/long grid. D show temperature, given again in Figure 4b? Can be avoided here.\nFigs. 2d and Fig. 4b: Gives the same information, only once for a specific RTS and once for all RTS. Maybe Fig. 2d is redundant (only mentioned in the discussion as an example where one could refer also to Fig. 4b)? At least have the same axis labels and legend.\nFig. 3: Nice figure, but not discussed, see comment above\nFig. 4: The figure is a bit messy, lots of different information, do you need the big pie diagrams? They more or less are similar. D is very messy, what do you want to show? High loss rates for arctic lowlands are never discussed really.\nFig. 5: This figure may be combined with Fig. 4b because of redundancies. The trend in 5b is really r>0.5? To be honest, for us the distribution of the red dots is quite random.\nFig. S2: We find the explanation of the figure is not clear: We assume that the red area around P1 is the deposition area? Why is the deposition earlier (2011) than the retreat of the headwall (2012-2015, indicated by Fig. S2g)? This is also corroborated by Fig. S2f, where deposition was in 2011 but head wall retreat until 2016? How can this be explained? Is it because only the largest elevation difference is analysed (see Line 1281 ff)? Then maybe mention it briefly here and not only in the supplement. Finally in Fig. S2b and e: Why are there so many elevations in 2015 available and only one for the other years?\nFig. S11a: The blue dots in a, are those all RTS? 2747? This is not possible.\nTable 1: Maybe add a column with loss rate. The loss rate is shown in Fig. 2a but it would be nice to have a summary for each region.\nOther comments:\nLine 37: Add the exact time period 2012-2022\nLine 98: You state that the existing literature focus on one or two-dimensional detection and measurement of RTSs, but give only an example for 2D detection. Please add another one for 1D detection.\nLine 162: Information in the brackets is unclear: Why give it in ha/ 1000 km2? This is not really intuitive. And what do the percentages describe? Many of such numbers in the text decrease readability, better to use a table to summarise numbers.\nLine 177ff: RTS and ice-cored moraines. We think there is quite a difference between recent ice-cored moraines like on Svalbard and old buried glaciers ice like on Banks Island (see major comment above). You identified 11 such RTS on Svalbard, this is quite few, considering the widespread distribution of ice-cored moraine complexes in Svalbard and certainly other Arctic islands.\nLine 222 ff: The rate is not really described in the results, this paragraphs mentions it but then goes straight to the acceleration. However, the rate could be interesting, see Fig. 4d: Why do the Arctic Lowlands have such a higher rate?\nLine 222 ff: How much sense does the analysis of mean acceleration make if the pattern show such pronounced peaks (see Fig. 4b)\nLines 321 and 336: the carbon loss is first given for decade, later again annually, confusing as only the exponents change when converting decadal to annual values. Suggest deleting the number in line 336.\nLine 378: What about uncertainties regarding the ERA5 data? Uncertainties regarding SOC as only upper 3m are used? Uncertainties in SOC database in general? (see also major comment above)\nLine 418: Is the mentioned uncertainty of 0.3-0.5m from literature or own validation?\nLine 420: Concerning the removed DEMS, did this have an impact on the overall analysis and results?\nLine 435: Why “e.g. larger than 3 times the RMSE”? Are there other significance criterions?\nLine 438: How is the uncertainty calculated?\nLine 444: What is b? Explanation not clear if reader is not familiar with Dai et al (2022).\nLine 463 ff: How were the numbers set for the classification of the confidence level?\nLine 478: Why criteria 3, is it necessary?\nLine 602: June 2012?\nLine 613 ff: Can the area-frequency distribution just be extrapolated to smaller RTS? Uncertainties? You can see that in Fig. S11a, the linear fit does not work for large areas. Why should it work for small ones?\nLine 745: Link for RTS shapefiles not found (see comment above)\nLine 1292: Do you mean “exclude winter scenes”?", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_4": "See also uploaded doc file, with same content.\n\nGeneral comment\nI am not at all a specialist of neither remote sensing nor deep learning methods. Therefore my comments will concentrate on general considerations and formal issues.\nThe paper offers an overview and pan arctic quantification of volume loss and soil organic carbon loss by retrogressive thaw slumps obtained through deep learning on a comprehensive DEM dataset.\nCompared to the already very rich literature on the subject, partly produced by the author team members themselves, the novelty of the contribution can be summarized as follows:\n- First evaluation of the phenomenon on a pan-arctic scale (previous works were done on regional levels)\n- First quantification of soil loss volumes on a decadal and annual time scales (previous works focussed on areas)\n- First estimation of soil organic carbon released by these slumps\n- The systematic use of the ArcticDEM data set.\nThe methodological approach seems to be robust. It uses a high resolution dataset, robust and recognized co-registration and extraction algorithms, several field evidence datasets for calibration and validation. The discussion on uncertainties, error sources and limitations is honestly done (text and mainly supplementary information) and gives a realistic idea of the reliability of the results.\nThe methodological developments (online methods) shows however that the approach needs a whole set of adjustments and model runs, and still gives approximative results. The results still need a huge work of visual/manual correction, which was done by the team. From this point of view, the paper should possibly be regarded rather as a methodological development paper than for its quantitative results.\n\nSome specific comments.\nTitle\nThe title doesn’t really reflect the content, and is written in too general terms. I suggested a more precise wording, like:\nPan-arctic quantification of total decadal and annual volume loss and soil organic carbon release by large retrogressive thaw-slumps during the decade 2012-2022.\nIn fact, the work doesn’t quantify the total mass loss, bur only the retreat during the considered decade. And it is explicitly restricted to retrogressive thaw slumps, excluding other mass wasting processes. This should be reflected in the title.\nArea estimates – lines 161-164\nSignificance of percentages given in brackets is unclear: … in areas with high ground ice content (1.37 ha 1000 km-2; Fig. S10; > 20 %), … etc.\n20 % of what ? ice-content ? of the data set ? … ?\nDiscussion – RTS Activity in response to Climate Trends – lines 252-280 + fig 4 and caption + fig 5 caption + fig S13\nText gives thresholds at 0.5, 0.3, -0.3, -0.5 for both temperature and precipitation. Same in figure 4 caption, and figure 5 caption\nBut in figure legends of figure 4, figure 5 and figure S13, the thresholds are given at 0.6, 0.2. -0.2, -0.6.\nCorrect either the figures or the text and caption, according to the thresholds used for the pie charts.\nSoil organic carbon release\nThe discussion on possible much lower SOC release than previously modelled is interesting. However, after reading the methodological part, and regarding the large limitations and uncertainties, it seems that your quantification rather underestimates the total volumes (missing landforms, limited considered depth, …) and should therefore be considered as a minimum (but I may have misunderstood the final result). May be this should be made more clear and less affirmative ?\nSources of error and analysis limitations\nThis sub-title is a bit peculiar at this place in the manuscript ! The following text is neither a conclusion nor a real uncertainty analysis.\nThis part should rather be placed right after the results, and before the discussion. It should best integrate the whole text S1, which is important for understanding.\nFigures\nMost figures are OK and informative.\nFigure 4\nThe figure 4a is small and contains too much information. Remove single dots and leave only the pie diagrams ?\nFigure 5\nFigure 5b: is the linear fit really significant ? Visually there is no fit at all …\nSupplementary material\nGood and useful for understanding.\nThe part on “Sources of error …” should be more explicitly referred to in the text, since it is essential for understanding the limitations of the study.\nReferences and refence list\nThere is unsufficient check of citations and reference list:\n• Many references cited in the text are missing in the reference list (highlighted red in the edited PDF).\n• For several references there is a wrong date, or an ambiguous date (e.g. missing precision for 2022a, 2022b) (year highlighted red in PDF).\n• Many references in the list are not cited. They have to be removed from the list (highlighted in red in the list).\n• Formatting of references in the list: many journal titles are not in italic (highlighted green in the PDF).\n• Some publications years do not correspond between text and reference list (year highlighted purple in PDF).\n• Some references near the end of the list are not in alphabetical order (highlighted with purple line number in PDF).\nPhilippe Schoeneich"}, "author_response": "The manuscript by Dai et al. is a highly ambitious piece of work that presents novel information on the dynamics of RTSs at the Pan-Arctic scale. By using large amounts of ArcticDEM data, the analysis covers the entire Pan-Arctic at a very high spatio-temporal resolution. The approach shares aspects (using ArcticDEM at Pan-Arctic scale with deep learning approaches) with the work by Huang et al. (2023 ISPRS Journal of Photogrammetry and Remote Sensing, https://doi.org/10.1016/j.isprsjprs.2023.10.008). However, the present manuscript goes beyond identifying RTSs by presenting spatially and temporally detailed volumetric estimates of soil loss and associated carbon release. The found magnitude of carbon release from RTSs provides an important new viewpoint to the globally relevant issue of permafrost degradation and associated carbon budgets. The manuscript also makes an interesting point by showing how volumetric RTS activity demonstrates latitudinal and regional differences in their responses to climate, but also that ecosystem properties markedly complicate this connection.\n\nA comment that needs to be made is regarding the relatively weak classification precision of the presented deep-learning method. High number of false positives is typical for approaches using automated recognition of landforms, and I appreciate the authors for addressing this challenge and pursuing to improve the automated mapping by combining solutions (change detection algorithm and deep learning) and using additional training data with data augmentation methods. Yet, also after this model precision remains low and coming up with reliable results requires undoubtedly laborious visual inspection.\nResponse: we carried out the evaluation in another test area in an Arctic desert environment with little vegetation, which gives a better precision (65%). Now we add the results to Table S3 and a paragraph to discuss this in the Methods section (Step 3b). We further note that the inventory was fully quality-checked by permafrost experts so that the samples included in the inventory are of the highest confidence.\nOn the other hand, the outcomes of the visual inspection are also affected by which locations are selected as candidate RTSs with the initial model and augmented models. Please see my specific comments regarding this aspect below.\nResponse: We now add: “Furthermore, the impact of the RTS missed in the deep learning detection is likely around 20% (recall is high, 77% to 80% - Table S3).”\nThe results are interesting and novel and should be of broad interest to the permafrost research community. The manuscript has an easy-to-follow structure even though the methodology presented involves several processing steps. The figures are clear and interpretable, and the language seems to me (as a non-native speaker) to be in very good order. Below I raise a couple issues and add minor suggestions for corrections and improvements. I recommend publication of the article after my raised issues have been addressed.\nSpecific comments:\nrr166-183: Here, the authors describe how the mapped RTS distributions are consistent with previous studies. I would also like to see a discussion of whether there are any regions or environments where RTS activity could not be mapped with a similar precision, or if the results contradict previous findings. Now, the discussion of uncertainties centers on the large number of false positives but the potentially missed RTS activity is less addressed. Although, the limitations in detecting some established activity at megaslumps due to narrow horizontal scale of activity was satisfactorily explained, which indicates that the approach may not have high recall in all conditions. Do we know if model evaluation statistics (from prediction and evaluation with Victoria-Peel Plateau data) hold for the Pan-Arctic area or can there be regional/environmental variation?\nResponse: we carried out additional evaluation in the Arctic desert environment, which gives a good precision of 65% and high recall of 77%.\nWe now add: “The prediction performance shows regional variations due to different land cover and environmental settings such as vegetation cover, snow and ice cover, and water bodies. Further evaluation is carried out in Eureka, Canada (Fig. S5c,56), which represents the Arctic desert environment with little vegetation. There are 45 RTS as ground truth in this test area manually examined by subject experts. Due to minimal impact from vegetation changes, the ensemble model achieves a considerably higher precision of 65% and a recall of 77% (Table S3).”\nTo address potentially missed RTS we add: “the impact of the RTS missed in the deep learning detection is likely around 20% (recall is high, 77% to 80% - Table S3).”.\nWe also discuss additional visual examples in the Discussion with locations from the literature of regions and environments where RTS activity could not be mapped because of our minimum size and elevation change threshold.\nr517: The split ratio of 9:1 sounds reasonable. Were other ratios tried in model training to examine their effect on overall accuracy?\nResponse: The split ratio 9:1 seems to be the widely used ratio. No other ratios were tested.\nrr563-564: Were the samples of 1660 training and 189 validation images also here split randomly based on the manually selected and inspected 1660+189 images? Please elaborate because models should be affected by how well the variation in the modelling data is represented.\nResponse: Yes, they are split randomly.\nNow we have: “The dataset is randomly split into training and validation sets at a 9:1 ratio, resulting in 1660 training images and 189 validation images (Fig. S7). ”\nrr576-578: I find it quite surprising how the combination of five models (original + augmented) increase precision so drastically. Model ensembles often do outperform single models, but the improvement is surprisingly large. Maybe single models identify different FPs and when put together, there are not so many common FPs, and thus the precision is higher? Could this be elaborated?\nResponse: We now add Fig. S6b to illustrate the improvement from ensemble modeling. The improvement (false positives from 357 to 210) is not very drastic, now we soften the tone in our sentences.\nNow we have: “since individual models may detect different false positives, when combined, the overlap of common false positives reduces, hence enhancing the precision (Fig. S6b). For the test area in Peel Plateau, the ensemble combination of five models (Table S3) reduces the number of false positives from 357 to 210.”.\nr745: I could not access the produced data via the (supposedly) DOI link (Mendeley?). Are the data available already?\nResponse: The data link was not ready. It should be ready now for download in the link provided in the data availability section. The associated metadata is still under review at Arctic Data Center.\nFig. S7-8. The final detected RTSs are seemingly concentrated in the training data (both 1 and 2) locations (as should be expected) while relatively little RTSs are detected in some areas known to have abundant occurrences. For example, only 13 RTSs have been detected in Alaska. I was wondering why there were no training images in Step 3b selected from Alaska at all? Could this be because the initial model from Step 3a could not predict potential RTSs there due to different environmental conditions in those areas, and consequently that also Step 3b results in very little RTSs in the western part of Brooks Range, for example? Or are these possibly “missing” RTSs simply filtered out in the classification in Step 2 for being too small or not active during the study period?\nResponse: No training images in Alaska were selected because the number of candidates in Alaska is too high, which is laborious for visual inspection. As shown in table S4, there are 19,154 candidates from step 3 (multimodel ensemble) for North Alaska. The number of candidates from step 3a (model 1) is normally higher, hence they are not missed in Step 3a.\nIndeed, a lot of slumps in the western part of Brooks Range are filtered out in Step 2. We checked and now visually address in the manuscript the published slumps in (Swanson and Nolan, 2018; Swanson, 2021). Two of those are detected in our study. The rest is either too small or filtered out in Step 2 because no pixels within the scarp area exhibited a significant elevation change exceeding 2 m, thereby failing to meet Criterion 2.\nWe added Fig. S14h to illustrate these cases. And we also added further clarification at the beginning of the results section and elsewhere throughout the document about our minimum sampling unit (minimum >10,000 m2 and minimum elevation change of 2 m) and make further direct statements such as “ Due to the minimum size or elevation change threshold adopted in this work, our dataset also excludes many small or shallow slumps characterized by smaller elevation changes (e.g., Fig. S15h;26,70).”\nTechnical touches:\nr99. Replace eg with e.g.\nResponse: corrected.\nFig. 1 a, b. Imagery source not visible, please enhance.\nResponse: We have added the following text at the end of the figure description to clarify the image sources: ESRI basemap imagery was used for (a), (b) and (d) accessed via ArcGIS Pro v2.7, Sources: Esri, DigitalGlobe, GeoEye, i-cubed, USDA FSA, USGS, AEX, Getmapping, Aerogrid, IGN, IGP, swisstopo, and the GIS User Community\nr223: 3gh -> 3g,h\nResponse: Here Fig. 3h is more relevant. We change it to Fig. 3h.\nTable 1: The pink and grey colors are hardly separable. Consider adding contrast.\nResponse: Changed as suggested.\nFig. 4: The legend for the correlation coefficients has different ranges to those presented in the figure caption and elsewhere in the text. Which are correct?\nResponse: The legend in the figure is now updated to be consistent with the caption.\nrr355-360: This sentence could be more clearly formulated. Is there a comma missing? Consider revising.\nResponse: We decide to retain this sentence.\nr643: Please elaborate the idea behind selecting the growth rate of double the long-term rate in Segal et al. 2016?\nResponse: Segal et al determined the mean standardized/normalized growth rate of 0.67% by considering changes between 1950-1960 imagery and 2004-2008. This growth rate is too low to reflect 21st century conditions as many studies from various areas have reported a recent intensification of thaw slumping. These intensifications are often reported in the literature as exponential, but to remain conservative in our assumptions we simply doubled Segal’s standardized growth rate to 1.3%.\nWe clarified the sentence by providing more details: “We multiplied the scar area estimates with a mean standardized growth rate of 1.34% (i.e., twice the long-term rate observed between 1950-1960 and 2004-2008 in 55 to account for RTS intensifications in the last two decades) to obtain an estimate of annual area change.”\nr652: “…for a pan-Arctic total of ±12,514”. Please check if this is correct. Are confidence intervals/standard deviation for the value missing here or how should one read the value?\nResponse: these should be read as one value. The plus/minus was removed in the text (and was supposed to indicate an approximation, not a confidence interval).\nr940: Please check the text font.\nResponse: Font corrected.\nFigs. S7: This is a bit rough figure. Consider using the same longitude of origin for the polar projection that was used in the other maps for better comparability. Please, also acknowledge the imagery source and consider adding a legend in which training and validation images are preferably presented with unique colors.\nResponse: Figure revised in polar stereographic projection. Google Earth image is acknowledged.\nReviewer #1 (Remarks on code availability):\nUnfortunately, I only found the Github link now as I am submitting my assessment. Anyway, I use different coding languages in my work, so I would not be able to test the reproducibility of the code or provide a sufficient assessment.\nResponse: Now we add a code availability section and instructions.\nReviewer #2 (Remarks to the Author):\nThe manuscript by Dai et al quantify retrogressive thaw slumps (RTS) in the northern permafrost area in terms of location, time of development and volume. The data is generated by comparing multitemporal DEMs based on ArcticDEM and identifying RTS based on a machine-learning algorithm. The data set not only identified polygons of RTS, but also the volume of the landform. This in turn could be used to estimate possible carbon release due to permafrost thaw and erosion. This combination makes the manuscript and its data base highly valuable for the research community, and this study is as far as we know unique and absolute worth publishing, even if the time series of 10 years is too short to evaluate real climate trends.\nWe are not experts on the machine-learning algorithms and let this part be evaluated by other reviewers. However, the following major comments are given and may be considered:\n• Title: Be precise, you deal (at least according to your writing) almost exclusively with retrogressive thaw slumps (RTS), why write then “large permafrost-thaw driven mass wasting events”? See also comment below.\nResponse: we now change the title as “Volumetric Quantifications and Dynamics of Areas Undergoing Retrogressive Thaw Slumping in the Northern Hemisphere”.\n• Type of RTS: We are not always sure what landforms you have detected. Of course, the big RTS in Arctic lowlands are captured and exemplified in images in the paper. However, in the Arctic islands, e.g. on Svalbard, you mentioned large RTS developed (by thermal erosion?) in large recent (< 250 years old) ice-cored moraine complexes in front of recent glacier fronts. Those those would have a quite different development history than ground-ice induced thaw slumps, even in areas of buried old Pleistocene glacier ice (there is a large literature base on this thermal erosion in ice-cored moraines on e.g. Svalbard). Maybe, you should show different types of identified RTS from the various regions? We wanted to have a look in your data for e.g. Svalbard, but the link provided failed.\nResponse: We now add a Fig. S13 to show examples of different RTSs such as thermoterraces and glacier dead-ice back-slumping in Svalbard. The data link should work now.\nWe also made further clarifications with respect to the scope of the study and the types of thaw-driven slumping we included.\n“Other areas undergoing RTS were detected along current glacier margins in northern Greenland, Svalbard, near Klutlan Glacier, Canada; and in Novaya Zemlya, Russia, where glacier dead-ice backwasting occurs (Fig. S13f) (e.g., 60). There are considerable geomorphological differences between dead-ice backwasting in recent ice-cored moraine complexes and RTS, yet the small number of large dead-ice backwasting detections (> 10,000 m2) were not manually excluded from the inventory as both mass wasting phenomena occur on ice-cored slopes, are thaw related, and exhibit similar elevation change “rainbow” patterns due to their retrogressive, rotational sliding failure mechanisms39–41. Furthermore, the strong linkage between RTS detections with modern-day glacier margins is compatible with RTS activity across northwestern Canada where terrain encompassing the former glacial margins of the Laurentide Ice Sheet contains relict glacier ice preserved by permafrost61.”\n• Correlation of events to temperature and precipitation: First of all, the results of the correlation analysis, which is interesting, is in the Discussion (chapter 3.1.), but it should be a Result. Consider moving the first part of that chapter (l. 252. onward).\nResponse: The first part Section (RTS activity in Response to Climate Trends) is now moved to Results.\nYou then write that 60% of the RTS are neither correlated with trends in ERA temperature or precipitation. I am not sure how you made this correlation, but an event one year is not necessarily correlated with climate the same year, but with the weather development the previous year. E.g. an extreme summer one year may lead to weakening of the structure and failure in one of the subsequent years. You should consider this in your discussion (se pre-conditioning papers in high-mountain permafrost slopes, e.g. here https://iopscience.iop.org/article/10.1088/1748-9326/ad3c55). This pre-conditioning may be more important further south than in the northernmost parts?\nResponse: We thank the reviewer for pointing this out, and in turn we carried out additional analyses concerning the correlation with the climate (temperature and precipitation) from the previous year. Indeed, the preconditioning effect is more prominent in the south than in the north regions. We added the results to a new Fig. 5a.\nWe also fully updated the results and discussion text considering same-year and previous year preconditioning effects.\n• Soil organics release: These are nice data, again, why do you not present these as a result, and then discuss them later.\nResponse: The section (Soil organic carbon release) is now moved to Results.\nLine 341 ff: The reference Tanski et al. (2017) does not support the statement of less SOC release. Here, the numbers are just small because SOC is buried again. However, in this study, the deposition of SOC is not taken into account, so that your numbers should be compared only to the released SOC.\nResponse: We have removed this reference from the manuscript. In places were field studies are used to support our argument that modelling studies have overestimated carbon release amounts we clarified our language used by specifying carbon released and net carbon amounts, modified text reads: “While our study shows less overall carbon release from RTS process, field studies19,20 on the Peel Plateau have shown that released carbon is sequestered within redistributing and accumulating sediments (e.g. debris tongues, Figure 1c) by becoming chemically bound through interactions with mineral soil. While variations in net carbon release vary regionally77, likely due to variation in carbon content of near-surface and deeper permafrost19, overall, our study along with these recent field studies19,20,77 point to net carbon released by RTS to be of lower magnitude than initially suggested by modelling18 and merits future work to better constraint values.”\n• Discussion in general: The Discussion in general could benefit from a more detailed elaboration of the results from Fig. 3, which is nice, but not discussed further in the paper (just presented as a result). E.g., why do you observe such areas, volumes, depths, rates and accelerations? Is air temperature and precipitation the only explanation as you focus on this in the discussion? What about ground ice contents, topography, geomorphology etc.?\nResponse: we added a new paragraph in the Discussion to elaborate on Fig 3 and what it means in the larger context of RTS studies. We also clarified our scope of this work as detection and description of patterns, but acknowledge that the drivers/reasons behind the areas/volumes/depths/rates/accelerations is an area of future work.\nSpecifically about the scope we added: “ It was beyond the scope of this study to determine the drivers behind the regional differences in volumes, depths, yields, scaling coefficients, or momentum, which should be considered in future work. Drivers governing these differences are most likely related to regional topographic slope configurations, subsurface conditions, landscape genesis (e.g., glacial legacy), and the most recent disturbance or thaw history (e.g., presence of old or ancient RTS disturbances, wildfire) (31,42,61,64,67). ”\n• Comment on the so-called “Rainbow pattern”: What about RTS where activity started or ended during the observation period and thus show no or only a weak rainbow pattern? Are these excluded? If that is the case it should be discussed in the uncertainties.\nResponse: We added a new paragraph in the Discussion to describe the uncertainties in more details, including weak rainbow patterns, RTS below the minimum size threshold or below the minimum elevation threshold. In the carbon release discussion we further added: “Nevertheless, the carbon release is likely underestimated because some RTS are missed due to a less prominent time transient change pattern if only a single year of retrogression occurred during the observation period or if the disturbance area was affected by surrounding noise.”\n• Uncertainties: A major result of the discussion is the relationship between RTS activity and climate, and the influence of RTS on carbon release. This discussion is based on ERA5 data and SOC map by Hugelius et al. Here, we miss a more in-depth description of the uncertainties in ERA5 and the SOC map, which of course influence the conclusions. soil organics carbon map.\nResponse: For the SOC results, we have added the following text to the discussion to clarify uncertainties:\n“Our calculated value of approximately 1.95×10-3 Pg carbon per year within our time period, which is the first attempt at this spatial scale with observed volumetrics, is much less (~0.2%) compared to what was estimated (using modeling) to be released annually from gradual thawing permafrost (~1 Pg carbon per year)11. It could be argued that a large number of small disturbance areas undergoing RTS (area < 10,000 m2, small elevation changes) not considered in this study could lead to underestimations of carbon release. To estimate this effect, we used area-frequency distribution modeling and comparisons with detailed manual RTS inventories to estimate the total RTS population size and volumetrics of small areas undergoing RTS (Figs. S11a, S12), resulting in an additional 2,244 to 12,514 areas or 18.0 to 63.8 ×106 m3 (see Methods). The inclusion of small RTS increases the total Northern Hemisphere volume estimate from 317.0 ×106 m3 to a range of 335 to 381 ×106 m3, and shows that the effect of undetected smaller RTS on carbon release is small, comprising only 5-17% of the estimated total volume. Furthermore, the impact of the RTS missed in the deep learning detection is likely around 20% (recall is high, 77% to 80% - Table S3). Nevertheless, the carbon release is likely underestimated because some RTS are missed due to a less prominent time transient change pattern if only a single year of retrogression occurred during the observation period or if the disturbance area was affected by surrounding noise. The mean depth of elevation loss is 3 m, so our estimates represent at least average conditions within the study domain. 58 found that the depths exceeding 3 m contributed to 15 to 20 % volume change and 5 to 10 % of total soil carbon mobilization. In our dataset, excluding carbon at depths greater than 3 m represents 30 % of volume with no carbon estimate. Ground ice content was not included in our carbon release estimate because the uppermost layers are comprised of the active layer and ice-poor permafrost above massive ice19. We believe the underestimation of not including the deeper ground layers offsets this potential overestimation for RTS with ground ice within the upper 3 m of soil.”\nFor uncertainties of ERA5, we have: “Limitations also arise from the heterogeneity in the accuracy and representativeness of ERA5 trends due to the limited number of Arctic weather stations and the short amount of yearly measurements (n ~10, Fig. 2d) used in this study. ”\n• Figures: Some of the figures are a bit messy and should be revised. A list of suggestions are given below in the minor comments.\nMinor comments are as follows:\nFigures and tables:\nFig. 1: Scales are lacking in a, c and d. Lat/long are missing in at least a and b.\nResponse: Scales are added to a, b and d. Lat/Long information has been added in the figure description: Examples of slump landforms in the Eureka Sound Lowlands (~79 °N) that encompass south-eastern Axel Heiberg Island and west-central Ellesmere Island, Canada….The example in (b) can be found at 79.49 °N; 86.84 °W and (d) at 79.61 °N; 86.23 °W.\nIs c and d the same RTS? If yes, which one of the RTS shown in b?\nResponse: No the field example is not the same as the remote sensing examples but they occur in the same area, this clarification has been added to the figure description: (c) shows a field example of a slump (different from those shown in (b) and (d), but from the same area) taken from a helicopter in 2016.\nFig. 2: Scale is missing in a, no lat/long grid. D show temperature, given again in Figure 4b? Can be avoided here.\nFigs. 2d and Fig. 4b: Gives the same information, only once for a specific RTS and once for all RTS. Maybe Fig. 2d is redundant (only mentioned in the discussion as an example where one could refer also to Fig. 4b)? At least have the same axis labels and legend.\nResponse: Scale and lat/long are now added.\nWe prefer to keep Fig. 2d, as it shows the correlation for individual RTS, as a nice transition to Fig. 4a.\nFig. 3: Nice figure, but not discussed, see comment above\nResponse: Fig. 3 is now discussed in the Discussion.\nFig. 4: The figure is a bit messy, lots of different information, do you need the big pie diagrams? They more or less are similar. D is very messy, what do you want to show? High loss rates for arctic lowlands are never discussed really.\nResponse: We deleted the pie diagrams from the figure. In the Results and Discussion we expanded on the loss rate information. For example, we provide explicit reasons as to why the Arctic Lowlands exhibited a higher volume loss rate (due to its large RTS population).\nFig. 5: This figure may be combined with Fig. 4b because of redundancies.\nResponse: We now add the correlation with the previous year temperature and precipitation to Fig. 5a, which shows the unique latitudinal variations in their responses to climate, showing different information from Fig. 4b.\nThe trend in 5b is really r>0.5? To be honest, for us the distribution of the red dots is quite random.\nResponse: We clarified in the discussion that the direct correlation is r = 0.23, and stated “a positive yet weak correlation was observed…”. We further added an error bar for the fit in Fig. 5b. “\nFig. S2: We find the explanation of the figure is not clear: We assume that the red area around P1 is the deposition area? Why is the deposition earlier (2011) than the retreat of the headwall (2012-2015, indicated by Fig. S2g)? This is also corroborated by Fig. S2f, where deposition was in 2011 but head wall retreat until 2016? How can this be explained? Is it because only the largest elevation difference is analysed (see Line 1281 ff)? Then maybe mention it briefly here and not only in the supplement.\nResponse: The earlies retreat in Fig. S2g started in 2011 (dark blue areas). Yes, the largest deposition seems to occur in 2011. The deposition from later years may flow and spread out downstream.\nIn the caption, we now add: “The deposition at this point (P1) appears to have occurred in 2011, while the deposition from later years may have flown and spread out downstream.”\n“The estimated time identifies when the most substantial elevation change occurred within the time series.”. We also added additional narrative about retaining only the largest elevation difference in the Discussion (second paragraph).\nFinally in Fig. S2b and e: Why are there so many elevations in 2015 available and only one for the other years?\nResponse: For local areas, our dataset (ArcticDEM) may have more data in some years and less data in other years. But generally, it has data coverage in most years. This is discussed in the text, e.g., we have: “an overall limitation of ArcticDEM is that the temporal resolution varies throughout its geographic coverage depending on the availability of stereographic pairs of images used to create the DEMs. Therefore, the temporal and geographic coverage of DEMs is not uniform across the study region.”.\nFig. S11a: The blue dots in a, are those all RTS? 2747? This is not possible.\nResponse: Yes, blue dots represent the histogram of 2747 RTS. Each dot is not an RTS, rather represents the probability density for RTS in each area bin (bin width 1000m2).\nNow we add: “𝑝 = 𝑐/(𝑁𝑤), where c is the number of RTS within an area bin, w is the bin width (1000 m2), N is the total number of RTS (2,747)”\nTable 1: Maybe add a column with loss rate. The loss rate is shown in Fig. 2a but it would be nice to have a summary for each region.\nResponse: A column of total volume loss rate is now added to Table 1.\nOther comments:\nLine 37: Add the exact time period 2012-2022\nResponse: Changed as suggested.\nLine 98: You state that the existing literature focus on one or two-dimensional detection and measurement of RTSs, but give only an example for 2D detection. Please add another one for 1D detection.\nResponse: we’ve provided additional references to support this statement.\nLine 162: Information in the brackets is unclear: Why give it in ha/ 1000 km2? This is not really intuitive. And what do the percentages describe? Many of such numbers in the text decrease readability, better to use a table to summarise numbers.\nResponse: we clarified what the percentages describe (ice content). There are different ways of providing density information, with RTS areas usually given in hectares and study areas in squared km. Also, because the total RTS area across the Arctic is small (and the study area is large), we felt that “1.37 ha 1000 km-2” was more intuitive than using the same area units, e.g., 0.0137 km2 1000 km-2.\nLine 177ff: RTS and ice-cored moraines. We think there is quite a difference between recent ice-cored moraines like on Svalbard and old buried glaciers ice like on Banks Island (see major comment above). You identified 11 such RTS on Svalbard, this is quite few, considering the widespread distribution of ice-cored moraine complexes in Svalbard and certainly other Arctic islands.\nResponse: Our identified number is small probably due to the area threshold (>10,000 m2).\nWe now add more notes in the Table S4 to clarify this: “ This table only includes medium- to large RTSs with areas larger than 10,000 m2.” and further expanded on the Results section:\n“Other areas undergoing RTS were detected along current glacier margins in northern Greenland, Svalbard, near Klutlan Glacier, Canada; and in Novaya Zemlya, Russia, where glacier dead-ice backwasting occurs (Fig. S13f) (e.g., 60). There are considerable geomorphological differences between dead-ice backwasting in recent ice-cored moraine complexes and RTS, yet the small number of large dead-ice backwasting detections (> 10,000 m2) were not manually excluded from the inventory as both mass wasting phenomena occur on ice-cored slopes, are thaw related, and exhibit similar elevation change “rainbow” patterns due to their retrogressive, rotational sliding failure mechanisms39–41. Furthermore, the strong linkage between RTS detections with modern-day glacier margins is compatible with RTS activity across northwestern Canada where terrain encompassing the former glacial margins of the Laurentide Ice Sheet contains relict glacier ice preserved by permafrost61.”\nLine 222 ff: The rate is not really described in the results, this paragraphs mentions it but then goes straight to the acceleration. However, the rate could be interesting, see Fig. 4d: Why do the Arctic Lowlands have such a higher rate?\nResponse: We expanded the section and clarified specifically why the Arctic Lowlands exhibited higher volume loss rates.\nLine 222 ff: How much sense does the analysis of mean acceleration make if the pattern show such pronounced peaks (see Fig. 4b)\nResponse: We acknowledge the comment and have revised the statement with further clarification and linkage to the next section (we do want to note that we used momentum, not mean momentum in this analysis.\n“Temporal metrics provided insights into the annual dynamics of RTS. The majority (75%) of areas undergoing RTS yielded <10,000 m3 annually (<20,000 m3 89.6%), with median yields of 5,013 m3 (Fig. 3g). Summarizing annual yields across regions (named after Provinces and Divisions as defined by 63), RTS located in the Arctic Lowlands (e.g., Banks Island), Mackenzie Mountains, and Brooks Range (incl. Peel Plateau) exhibited the largest yields of 10.0, 4.1, and 4.1 million m3/year, respectively, due to their large population of areas undergoing RTS (Table 1). At a combined 18.2 million m3/year, these three regions contributed up to 56.8% of the observed annual yields across the Northern Hemisphere. The majority of RTS exhibited a relatively stable momentum in their annual yields (i.e., near 0 m3/year2), meaning that annual activity is neither accelerating or decelerating. However, as a whole Northern Hemisphere population, momentum was slightly positive with a mean of 179 m3/year2 (Fig. 3h) indicating a slight acceleration in annual volume loss. These momentum estimates are affected by pronounced annual peaks and lows of activity, but provide useful insights into the net long-term trends of RTS in the Northern Hemisphere as well as sub-regions.”\nLines 321 and 336: the carbon loss is first given for decade, later again annually, confusing as only the exponents change when converting decadal to annual values. Suggest deleting the number in line 336.\nResponse: We delete the number. And we now only use Pg carbon per year to be consistent.\nLine 378: What about uncertainties regarding the ERA5 data? Uncertainties regarding SOC as only upper 3m are used? Uncertainties in SOC database in general? (see also major comment above)\nResponse: SOC uncertainties and ERA5 uncertainties are addressed in the previous response.\nLine 418: Is the mentioned uncertainty of 0.3-0.5m from literature or own validation?\nResponse: Mainly from literature and our group’s other related work.\nLine 420: Concerning the removed DEMS, did this have an impact on the overall analysis and results?\nResponse: Deleting bad quality DEMs has little impact on the overall analysis since bad DEMs only occur occasionally.\nLine 435: Why “e.g. larger than 3 times the RMSE”? Are there other significance criterions?\nResponse: We listed five conditions for identifying confidence levels (significance) in the last paragraph of Step 1.\nLine 438: How is the uncertainty calculated?\nResponse: We now add: “Its uncertainty is estimated based on the posteriori variance of unit weight and error propagation (Dai et al., 2022).”\nLine 444: What is b? Explanation not clear if reader is not familiar with Dai et al (2022).\nResponse: Now we add: “a+b represents the surface elevation after the event.”\nLine 463 ff: How were the numbers set for the classification of the confidence level?\nResponse: The thresholds are empirical.\nLine 478: Why criteria 3, is it necessary?\nResponse: The criterion 3 is chosen to exclude cases where there is one single event within the scarp area. This criterion will help select slumps with multiple thawing events which exhibit the unique rainbow pattern. We now add: “Criterion 3 is designed to identify slumps with multiple thawing events that display the distinctive rainbow pattern.”\nLine 602: June 2012?\nResponse: It’s now changed to: “summer defined as the period from June 20th to September 22nd each year”\nLine 613 ff: Can the area-frequency distribution just be extrapolated to smaller RTS? Uncertainties? You can see that in Fig. S11a, the linear fit does not work for large areas. Why should it work for small ones?\nResponse: The area-frequency distribution function is commonly used in landslide studies to estimate the number of small, undetected landslides, so indeed it can be extrapolated for this first-order estimate. This is one estimation method we used, which yields the volume of small RTS is about 5-8% of the entire population. This estimation generally agrees with the second method in the next paragraph, which shows that the volume of small RTS is about 17% of the entire population.\nLine 745: Link for RTS shapefiles not found (see comment above)\nResponse: New link is provided.\nLine 1292: Do you mean “exclude winter scenes”?\nResponse: We include winter scenes. To clarify this, we now add: “we had to include winter scenes to ensure adequate temporal data coverage. ”\nReviewer #3 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nReviewer #4 (Remarks to the Author):\nSee also uploaded doc file, with same content.\nGeneral comment\nI am not at all a specialist of neither remote sensing nor deep learning methods. Therefore my comments will concentrate on general considerations and formal issues.\nThe paper offers an overview and pan arctic quantification of volume loss and soil organic carbon loss by retrogressive thaw slumps obtained through deep learning on a comprehensive DEM dataset.\nCompared to the already very rich literature on the subject, partly produced by the author team members themselves, the novelty of the contribution can be summarized as follows:\n- First evaluation of the phenomenon on a pan-arctic scale (previous works were done on regional levels)\n- First quantification of soil loss volumes on a decadal and annual time scales (previous works focussed on areas)\n- First estimation of soil organic carbon released by these slumps\n- The systematic use of the ArcticDEM data set.\nThe methodological approach seems to be robust. It uses a high resolution dataset, robust and recognized co-registration and extraction algorithms, several field evidence datasets for calibration and validation. The discussion on uncertainties, error sources and limitations is honestly done (text and mainly supplementary information) and gives a realistic idea of the reliability of the results.\nThe methodological developments (online methods) shows however that the approach needs a whole set of adjustments and model runs, and still gives approximative results. The results still need a huge work of visual/manual correction, which was done by the team. From this point of view, the paper should possibly be regarded rather as a methodological development paper than for its quantitative results.\nResponse: We understand the reviewers’ point-of-view but disagree with regarding this manuscript to be solely a methodological development paper. Most deep learning papers have results typically limited to describing accuracy and runtime. In this manuscript we document for the first time a representative quantification of volumetric change from large areas undergoing retrogressive thaw slumping in the Northern Hemisphere, whereby we contribute new insights into the magnitude and manifestation of hillslope thermokarst by means of volumetric sediment budgets, non-linear scaling, and regional patterns We provide context and interpretations based on transparent estimates and implemented considerable due diligence beyond what a typical remote sensing/methodological paper would provide.\nSome specific comments.\nTitle\nThe title doesn’t really reflect the content, and is written in too general terms. I suggested a more precise wording, like:\nPan-arctic quantification of total decadal and annual volume loss and soil organic carbon release by large retrogressive thaw-slumps during the decade 2012-2022.\nIn fact, the work doesn’t quantify the total mass loss, bur only the retreat during the considered decade. And it is explicitly restricted to retrogressive thaw slumps, excluding other mass wasting processes. This should be reflected in the title.\nResponse: The title is now changed to: “Volumetric Quantifications and Dynamics of Areas Undergoing Retrogressive Thaw Slumping in the Northern Hemisphere”\nArea estimates – lines 161-164\nSignificance of percentages given in brackets is unclear: … in areas with high ground ice content (1.37 ha 1000 km-2; Fig. S10; > 20 %), … etc.\n20 % of what ? ice-content ? of the data set ? … ?\nResponse: Now it is changed to: ice content > 20%; ice content between 10% and 20 %; ice content < 10%.\nDiscussion – RTS Activity in response to Climate Trends – lines 252-280 + fig 4 and caption + fig 5 caption + fig S13\nText gives thresholds at 0.5, 0.3, -0.3, -0.5 for both temperature and precipitation. Same in figure 4 caption, and figure 5 caption\nBut in figure legends of figure 4, figure 5 and figure S13, the thresholds are given at 0.6, 0.2. -0.2, -0.6.\nCorrect either the figures or the text and caption, according to the thresholds used for the pie charts.\nResponse: Legends in Figures 4 and S13 are now updated to be consistent with the text and caption.\nSoil organic carbon release\nThe discussion on possible much lower SOC release than previously modelled is interesting. However, after reading the methodological part, and regarding the large limitations and uncertainties, it seems that your quantification rather underestimates the total volumes (missing landforms, limited considered depth, …) and should therefore be considered as a minimum (but I may have misunderstood the final result). May be this should be made more clear and less affirmative ?\nResponse: acknowledged. We revised the discussion on this by further integrating the estimates of small slumps that felt outside the scope of our detection models.\n“It could be argued that a large number of small areas undergoing RTS (area < 10,000 m2, small elevation changes) not considered in this study could lead to underestimations of carbon release. To estimate this effect we used area-frequency distribution modeling and comparisons with detailed manual RTS inventories to model the total RTS population size and volumetrics of small areas undergoing RTS (Figs. S11a, S12), resulting in an additional 2,244 to 12,514 areas or 18.0 to 63.8 ×106 m3 (see Methods). The inclusion of small RTS increases the total Northern Hemisphere volume estimate from 317.0 ×106 m3 to a range of 335 to 381 ×106 m3, and shows how the predicted impact of undetected smaller RTS on carbon release is small (only 5-17% of the estimated total volume). Furthermore, the impact of missed RTS from the deep learning detection is likely low, around 20%, since the recall is high, 77% to 80% (Table S3).\nNevertheless, the carbon release might be underestimated considering that some RTS might be missed due to a less prominent rainbow pattern if it had only a single year of retrogression during the observation period or if it was affected by surrounding noises”.\nSources of error and analysis limitations\nThis sub-title is a bit peculiar at this place in the manuscript ! The following text is neither a conclusion nor a real uncertainty analysis.\nThis part should rather be placed right after the results, and before the discussion. It should best integrate the whole text S1, which is important for understanding.\nResponse: we now integrate and further expand on the whole text S1 in the discussion.\nFigures\nMost figures are OK and informative.\nFigure 4\nThe figure 4a is small and contains too much information. Remove single dots and leave only the pie diagrams ?\nResponse: We now remove the pie diagrams.\nFigure 5\nFigure 5b: is the linear fit really significant ? Visually there is no fit at all …\nResponse: We now add gray shading to show the standard deviation of the fit. In our discussion, we acknowledge the weak correlation, e.g, we have: “a positive yet weak correlation was observed…”\nSupplementary material\nGood and useful for understanding.\nThe part on “Sources of error …” should be more explicitly referred to in the text, since it is essential for understanding the limitations of the study.\nResponse: We now integrate the whole text S1 to the main text (Discussion).\nReferences and refence list\nThere is unsufficient check of citations and reference list:\n• Many references cited in the text are missing in the reference list (highlighted red in the edited PDF).\n• For several references there is a wrong date, or an ambiguous date (e.g. missing precision for 2022a, 2022b) (year highlighted red in PDF).\nResponse: Revised as suggested. See the following point-to-point response (For detailed comments in the PDF file).\n• Many references in the list are not cited. They have to be removed from the list (highlighted in red in the list).\nResponse: Revised as suggested\n• Formatting of references in the list: many journal titles are not in italic (highlighted green in the PDF).\nResponse: References formatting is revised to be consistent with the journal style.\n• Some publications years do not correspond between text and reference list (year highlighted purple in PDF).\n• Some references near the end of the list are not in alphabetical order (highlighted with purple line number in PDF).\nResponse: Revised as suggested\nFor detailed comments in the PDF file:\nLine 42: ‘Moreover’ -> ‘However’\nResponse: Revised as suggested.\nLine 71: Smith et al., 2023\nResponse: Reference added.\nLine 74, 83: IPCC,2023\nResponse: Citation corrected.\nLine 77: Ward Jones et al., 2022; Kokelj et al., 2023\nResponse: Corrected.\nLine 79, 90: Turetsky et al., 2020\nResponse: Reference added.\nLine 88: Sahkil et al., 2021\nResponse: Reference removed.\nLine 96: Zwieback et al., 2018\nResponse: Reference added.\nLine 97: Bernhard et al., 2022 a or b?\nResponse: Corrected as: Bernhard et al., 2022a\nLine 110, 177: Nitze et al., accepted\nResponse: Updated as Nitze et al., 2024.\nLine 115: Burn and Lewkowicz, 1990\nResponse: Reference added.\nLine 179: Kjær and Krüger, 2001\nResponse: Reference added.\nLine 263: Edlund and Garneau, 2000\nResponse: Reference added.\nLines 331-332: Turetsky et al., 2020\nResponse: Reference added.\nLine 343, 355: Shakil et al., (2020)\nResponse: Reference added.\nLine 395: Dai and Howat, 2018\nResponse: Reference is now added.\nLine 610, 1384: Welch, 1967\nResponse: Reference is now added.\nLine 634, 687: Bernhard et al., (2022) “a or b?”\nResponse: Corrected as: Bernhard et al., (2022a)\nLine 644: Segal et al.; (2016) “a, b or c?”\nResponse: It is Segal et al., 2016a. Corrected in the new citation.\nLine 683: Larsen et al., 2010; Tseng et al., 2013\nResponse: Reference is now added.\nLine 726: Wessel and Smith, 1991\nResponse: Reference is included.\nPhilippe Schoeneich"}, {"round": 2, "reviewer_comments": {"reviewer_1": "I wish to thank the authors for their careful consideration of my and other reviewers’ comments. I see that the manuscript is now improved in its methodological clarity and overall structure of the presentation. As a minor comment, I think that it could be considered if the text on the contribution of small RTS disturbance (<10,000 m2) to carbon release in the discussion should be already presented in the results section (“Soil organic carbon release”). Moreover, around lines 418-427 in the discussion there is a lot of repeating information from the Results, section “RTS Activity in Response to Climate Trends”. Apart from these notions, I recommend publication.", "reviewer_2": "Revision Review “Volumetric quantifications and dynamics …..”\nThe manuscript by Dai et al quantify retrogressive thaw slumps (RTS) in the northern permafrost area in terms of location, time of development and volume. The data is generated by comparing multitemporal DEMs based on ArcticDEM and identifying RTS based on a machine-learning algorithm. The data set not only identified polygons of RTS, but also the volume of the landform. This in turn could be used to estimate possible carbon release due to permafrost thaw and erosion. This combination makes the manuscript and its data base highly valuable for the research community, and this study is as far as we know unique and absolute worth publishing.\nWe now have received a revised version of the manuscript, and we are pleased to see that many of our suggestions from the initial review has been followed, and that most of our comments are incorporated. The paper appears now much better, with a new title, better division of results and discussion, a improved discussion of landforms, incorporation of a statistical analysis of previous years temperatures and major improvements of many figures. However, some issues remained:\n1. Download of the shape files: It is now possible that the individual RTS can be loaded into a GIS, but they are simply numbered, so you have to download each one individually to find the one for a certain region. Maybe the authors could divide it into subfolders with the individual regions. In addition, we couldn't load the shapefiles with the coordinates for all slides. Could possibly be due to the fact that e.g. the PRJ file could not be downloaded (opens only in a separate table).\n2. The discussion is lengthy, and readability would be improved if divided into subsections with own headings, e.g. “Limitations and uncertainties”, “Volumetric yield estimates”, “Relation to climate trends”.\n3. Concerning climate relationships, avoid the term “climatic trend” for a 10-year time series. Determining real trends in climate demands at least 30 years, and we are still critical to the significance of figure 5b, even if the figure has improved. It is also normal to give R2-values (not only simple r-values). R2 would give below 25% of explanation, which is not much. It does not mean that climate is not influencing, but you cannot show it statistically significantly on a certain level in your study.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_4": "The main issues raised by my review were ceottectly adressed in the new manuscript version.\nThis can thus be published.", "reviewer_5": "My comments focus on machine learning section only. Overall, the methodology presented in the manuscript is rigorous, reliable, and represents a well-justified technical approach.\n1. In this manuscript, deep learning primarily serves as a tool to reduce manual workload, utilizing existing adapted methods. The main novelty lies in integrating change detection and classification into a unified framework, with prior knowledge incorporated via three selected feature bands.\n2. For machine learning-based land surface information extraction, training sample distribution is critical. As noted in the manuscript, performance varies across land cover types, indicating limited spatial transferability of the method. Thus, training samples should comprehensively located in the study area.\n3. The two-stage iterative training aligns with active learning principles, which is meaningful. While model ensemble yields marginal improvement, the method (e.g., voting or other techniques) should be clarified. What is the size of sample blocks?\n4. Using 10% of samples for accuracy validation may slightly overestimate performance compared to ground truth, but this does not undermine the conclusions.\n5. The core identification relies on predefined rules (e.g., elevation change, area thresholds). Machine learning here functions as a labor-saving tool rather than a solution driver, so high methodological novelty is unnecessary."}, "author_response": "Reviewer #1 (Remarks to the Author):\nI wish to thank the authors for their careful consideration of my and other reviewers’ comments. I see that the manuscript is now improved in its methodological clarity and overall structure of the presentation. As a minor comment, I think that it could be considered if the text on the contribution of small RTS disturbance (<10,000 m2) to carbon release in the discussion should be already presented in the results section (“Soil organic carbon release”). Moreover, around lines 418-427 in the discussion there is a lot of repeating information from the Results, section “RTS Activity in Response to Climate Trends”. Apart from these notions, I recommend publication.\nResponse: Thank you for your constructive suggestions.\nWe now move the text about small RTS to the results section.\nAs suggested, we delete the repeating details: “…(from 40% to 88% of RTS binned by latitude with r >=0.5 as latitudes increased from 71°N to 82°N). Precipitation-driven RTS decreased linearly from 60% to 13% with increasing latitude from 60°N to 82°N.”\nWe keep a brief summary of results as a transition to the discussion of the explanation of the phenomenon.\nReviewer #1 (Remarks on code availability):\nI am not sufficiently qualified to assess the code in the provided format and language.\nReviewer #2 (Remarks to the Author):\nRevision Review “Volumetric quantifications and dynamics …..”\nThe manuscript by Dai et al quantify retrogressive thaw slumps (RTS) in the northern permafrost area in terms of location, time of development and volume. The data is generated by comparing multitemporal DEMs based on ArcticDEM and identifying RTS based on a machine-learning algorithm. The data set not only identified polygons of RTS, but also the volume of the landform. This in turn could be used to estimate possible carbon release due to permafrost thaw and erosion. This combination makes the manuscript and its data base highly valuable for the research community, and this study is as far as we know unique and absolute worth publishing.\nWe now have received a revised version of the manuscript, and we are pleased to see that many of our suggestions from the initial review has been followed, and that most of our comments are incorporated. The paper appears now much better, with a new title, better division of results and discussion, a improved discussion of landforms, incorporation of a statistical analysis of previous years temperatures and major improvements of many figures. However, some issues remained:\n1. Download of the shape files: It is now possible that the individual RTS can be loaded into a GIS, but they are simply numbered, so you have to download each one individually to find the one for a certain region. Maybe the authors could divide it into subfolders with the individual regions. In addition, we couldn't load the shapefiles with the coordinates for all slides. Could possibly be due to the fact that e.g. the PRJ file could not be downloaded (opens only in a separate table).\nResponse:\nWe now group individual RTS based on regions.\nWe also add a README.txt file, and add the ftp link for each RTS in ArcticRTS.shp for easy download, as shown in the screenshot below:\nTo download VolLatlon.prj file from https://arcticdata.io/data/10.18739/A22805140/RTS/ , right click VolLatlon.prj, and then click ‘Save link as…’\n2. The discussion is lengthy, and readability would be improved if divided into subsections with own headings, e.g. “Limitations and uncertainties”, “Volumetric yield estimates”, “Relation to climate trends”.\nResponse: Subsection headings have now been added.\n3. Concerning climate relationships, avoid the term “climatic trend” for a 10-year time series. Determining real trends in climate demands at least 30 years, and we are still critical to the significance of figure 5b, even if the figure has improved. It is also normal to give R2-values (not only simple r-values). R2 would give below 25% of explanation, which is not much. It does not mean that climate is not influencing, but you cannot show it statistically significantly on a certain level in your study.\nResponse: We change ‘climatic trend’ to ‘climate variability’.\nWe agree that the significance of Figure 5b is limited, and we now delete it and the related discussions.\nReviewer #3 (Remarks to the Author):\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nResponse: Thank you for your constructive suggestions.\nReviewer #4 (Remarks to the Author):\nThe main issues raised by my review were ceottectly adressed in the new manuscript version.\nThis can thus be published.\nResponse: We appreciate your recommendation.\nReviewer #5 (Remarks to the Author):\nMy comments focus on machine learning section only. Overall, the methodology presented in the manuscript is rigorous, reliable, and represents a well-justified technical approach.\n1. In this manuscript, deep learning primarily serves as a tool to reduce manual workload, utilizing existing adapted methods. The main novelty lies in integrating change detection and classification into a unified framework, with prior knowledge incorporated via three selected feature bands.\nResponse: Thank you for your positive feedback.\n2. For machine learning-based land surface information extraction, training sample distribution is critical. As noted in the manuscript, performance varies across land cover types, indicating limited spatial transferability of the method. Thus, training samples should comprehensively located in the study area.\nResponse: Thanks for pointing out the importance of training sample distribution. As shown in our Figure S7, our training data distribution indeed covers various reigons in the Arctic. We now add your comments to the Methods section: “The performance differences across land cover types indicate limited spatial transferability of our method.”.\n3. The two-stage iterative training aligns with active learning principles, which is meaningful. While model ensemble yields marginal improvement, the method (e.g., voting or other techniques) should be clarified. What is the size of sample blocks?\nResponse: Thank you for the comment. We now clarify the ensemble method by adding: “The ensemble model is retrieved by applying a logical AND operator to all individual model outputs.”.\nWe also add description about the size of sample blocks in the previous section (Step 3a) when we describe data preparation: “The image sample size ranges from approximately 10,000 m² to 500,000 m².”\n4. Using 10% of samples for accuracy validation may slightly overestimate performance compared to ground truth, but this does not undermine the conclusions.\nResponse: Thank you. We will consider exploring different validation sample size percentages in future work.\n5. The core identification relies on predefined rules (e.g., elevation change, area thresholds). Machine learning here functions as a labor-saving tool rather than a solution driver, so high methodological novelty is unnecessary.\nResponse: We agree with your comments.\nReviewer #5 (Remarks on code availability):\nThe authors provided a README document and a tutorial for reproducibility. The workflow primarily relies on cited third-party tools, with minimal custom-developed components. Consequently, I did not conduct in-depth testing or evaluation of the code.\nResponse: We appreciate your review of our code and your comments."}]

Decision letter