Detailed Characterization and Monitoring of a Retrogressive Thaw Slump from Remotely Piloted Aircraft Systems and Identifying Associated Influence on Carbon and Nitrogen Export

Abstract

Ice-rich permafrost landscapes are sensitive to ongoing changes in climate. Permafrost retrogressive thaw slumps (RTSs) represent one of the more abrupt and prolonged disturbances, which occur along Arctic river and lake shorelines. These features impact local travel and infrastructure, and there are many questions regarding associated impacts on biogeochemical cycling. Predicting the duration and magnitude of impacts requires that we enhance our knowledge of RTS geomorphological drivers and rates of change. Here we demonstrate the utility of remotely piloted aircraft systems (RPAS) for documenting the volumetric change, associated drivers and potential impacts of the largest active RTS along the Old Crow River in Old Crow Flats, Yukon, Canada. RPAS surveys revealed that 29,174 m³ of sediment was exported during the initial evacuation in June 2016 and an additional 18,845 m³ continued to be exported until June 2019. More sediment export occurred during the warmer 2017 summer that experienced less cumulative rainfall than summer 2018. However, several rain events during 2017 were of higher intensity than during 2018. Overall mean soil organic carbon (SOC) and total nitrogen (TN) within sampled thaw slump sediment was 1.36% and 0.11%, respectively. A combination of multispectral, thermal and irradiance (derived from the RPAS digital surface model) data provided detailed classification of thaw slump floor terrain types including raised dry clay lobes, shaded and relatively stable, and low-lying evacuation-prone sediments. Notably, the path of evacuation-prone sediments extended to a series of ice wedges in the northern headwall, where total irradiance was highest. Using thaw slump floor mean SOC and TN values in conjunction with sediment bulk density and thaw slump fill volume, we estimated that 713 t SOC and 58 t TN were exported to the Old Crow River during the three-year study. Findings showcase the utility of high-resolution RPAS datasets for refining our knowledge of thaw slump geomorphology and associated impacts.

1. Introduction

Ice-rich permafrost landscapes are sensitive to increases in temperature and precipitation. Retrogressive thaw slumps represent a particularly dramatic landscape response that is expected to intensify in magnitude and frequency with future climate change [1,2,3,4]. These thermokarst slope disturbances are highly dynamic and contribute large volumes of materials downslope to lakes, drainage networks and coastal zones [5,6]. Initiation of thaw slumps can be triggered by a combination of local meteorological and geomorphological conditions. For example, headward fluvial erosion along riverbanks or coastal erosion in maritime environments can initiate terrestrial detachment failures, especially during extended warm and wet conditions [5]. Exposure of ice-rich permafrost enables thaw slump development, which can persist for many years depending on meteorological conditions. The general geomorphic setting of an active thaw slump following initiation includes a near-vertical headwall, which contributes thawed or thawing material to the headscarp, where ablation continues, transporting moisture and sediment to the slump floor [7]. When conditions are sufficiently wet, thawed material can be transported by rill erosion, fluvial transport and mass flows [3,8,9,10]. Flow of water and sediment, during either fluvial or more deep-seated mud-flow events [3], results in accumulation of debris and chemical constituents at the debris tongue, where they are available for mobilization and deposition in lakes or along drainage networks.

There is a strong consensus that the size, frequency and activity of thaw slumps generally increase with increasing temperature and precipitation. However, thaw slump development can be highly variable among and within landscapes depending on local ground ice content, physiography and geomorphic conditions. The largest thaw slumps in north America, referred to as mega-slumps, have been observed in the Richardson Mountains and Peel Plateau regions, NT, Canada [11]. Thaw slump activity here and close by in the Mackenzie Delta has increased in response to wet summer conditions [1,3]. Wet and warm conditions have also been linked with recent increases frequency of thaw slump formation in the Qinghai–Tibet Plateau [12]. Observed thaw slump activity along coastal regions including Tuktoyaktuk Peninsula, Canada and the Bykovsky Peninsula, Russia have been attributed to the complex seasonal interplay of available energy and late-season rainfall [13]. Other Canadian coastal thaw slumps on Herschel Island have experienced accelerated activity during recent years characterized by warm temperatures [14]. In addition to meteorological conditions, interactions with marine environments, including thermo-abrasion from waves and ice, can have a strong influence on initiation and perpetuation of thaw slump activity along coastlands [15]. Studies within inland environments have highlighted the role that local geomorphological conditions have on initiation of detachment failures (e.g., river headward erosion), although meteorological conditions have been identified as the main drivers of sustained development and activity.

Orientation relative to available solar energy also represents an important characteristic for perpetuating thaw slump development. Lewkowicz [7] noted that the direction of thaw slump headwall retreat can be greater along south and southwest-facing slopes with relatively high incident radiation. Though the Batagay mega slump in the northern Yakutia region, Russia drains to the north, during recent decades it has experienced the greatest expansion along the south and southeast-facing headwall [16]. While thaw slumps experiencing northward headwall retreat have been highlighted in other studies in NT, Canada e.g., [10,17], orientation was identified as less of a control on headwall retreat in the Noatak Valley, Alaska [18]. It can be expected that local geomorphological conditions determine whether orientation of thaw slumps, relative to the maximum solar angle, is an important determining factor of the extent and direction of headwall retreat. Few studies have documented the growth rate directions of active thaw slumps relative to solar angles. However, this insight may be key for modelling thaw slump development and associated hazards in many permafrost environments.

Thaw slumps can have a strong influence on biogeochemical cycling in lakes [19,20,21] and rivers [22,23,24,25]. The biogeochemical effects of permafrost thaw on aquatic systems is highly variable and complicated by relief, ice content, permafrost extent, parent material [25] and flow pathway hydrological connectivity [26]. The thawed sediment and solutes mobilized to aquatic environments are typically rich in nutrients and ions. Carbon concentrations are particularly high in Yedoma permafrost and are available for efficient decomposition upon thaw [27,28]. As microbial activity ensues within thawed permafrost, carbon can be released as CO₂ or CH₄ [29] or be mobilized to adjacent and downstream aquatic environments. Partitioning of these pathways is required for effective modelling of associated climate-change feedbacks. However, resolving this necessitates insight of the interactions among permafrost thaw and local conditions (e.g., surface hydrology [27]). However, local geomorphic conditions may be generalized depending on the scope of research. For example, while studies have identified that soil organic carbon (SOC) concentrations can increase with depth depending on location [30] and according to landscape feature type [31], nutrient availability for mobilization from thaw disturbance are often assumed uniform at the landscape scale. There is limited insight of how varying geomorphic and physiographic characteristics within and among thaw slumps may be spatially associated with thawed soil nutrient properties, and biogeochemical pathways (e.g., mineralization and/or export of SOC). Developing refined partitioning of carbon mobilization pathways requires an integrated approach [32].

Here we couple high-resolution mapping products with in-situ measurements of sediment carbon and nitrogen to track thaw slump volumetric change during a three-year period and identify key geomorphologic characteristics that influence SOC and total nitrogen (TN) availability. The mapping products were compiled from images acquired using a remotely piloted aircraft system (RPAS). Note that we prefer using the term RPAS since it is gender neutral and implies the requirement of human expertise, unlike other terms that are often used interchangeably including unmanned aerial/aircraft systems/vehicles (UAS/UAVs). RPAS is also used internationally (e.g., by the International Civil Aviation Organization) and by Transport Canada. Using RPAS, we track thaw slump development, including export volume and headwall retreat rate, in response to varying interannual meteorology, and positioning and orientation of thaw slump geomorphic components relative to solar angles (i.e., sun altitude and azimuth). We utilize high-resolution mapping products derived from structure-from-motion (SfM) photogrammetric analysis of RPAS-derived photographs to distinguish the physiographic and geomorphic components of a thaw slump. In addition, we estimate the amount of SOC and TN export to the Old Crow River system based on coupled in situ measurements and monitored volume export. The contributions of this study include demonstration of leading-edge methodologies for tracking thaw slump geomorphic change and drivers, and identification of the influence that these processes have on soil nutrient properties and export to the drainage network.

Study Site

This study was focused on the largest and most active retrogressive thaw slump in Old Crow Flats (OCF), northern Yukon during the period 2016–2019. OCF is the traditional territory of the Vuntut Gwitchin First Nation and is internationally recognized for its ecological and cultural integrity by the Ramsar Convention on Wetlands (1982). To our knowledge, no research has been conducted on the dynamics of retrogressive thaw slumps in OCF. However, this important landscape is in relatively close proximity to other regions that have experienced an increase in frequency and magnitude in thaw slump activity during recent decades e.g., [4,14].

The retrogressive thaw slump of interest (hereafter referred to as Slump1) is located 38 km north of the town of Old Crow (85 km by boat) on the Old Crow River (Figure 1). With the surrounding area reaching elevation over 300 m (above sea level), Slump1 likely exposed and mobilized a combination of glaciolacustrine (Holocene) and Pleistocene sediments [33]. Like most other thaw slumps in OCF, Slump1 is located on the western bank of the Old Crow River and faces southeast. The initial detachment failure occurred during June 2016, approximately one month prior to our initial survey (29 July 2016). During the first summer, the debris tongue extending from Slump1 blocked approximately 60% (~30 m) of the river width.

During this study, meteorological conditions were recorded by the Government of Canada (https://climate.weather.gc.ca/) at the Old Crow Airport (Figure 2a,b). June 2016 experienced the greatest cumulative daily precipitation (31.4 mm) and frequency (n = 11) for all spring months during this study. This corresponds with the timing of the initial evacuation of Slump1. July 2016 had the second highest frequency (n = 14) and cumulative daily rainfall (48.6 mm). The following 2017 spring-summer seasons began relatively dry, with only one rain event recorded in June (8.3 mm) followed by seven events during July (28.7 mm). August 2017 experienced 13 rain events, which totaled 76.6 mm. June 2018 was wetter than June 2017, with 11.9 mm rainfall occurring over 13 events. As well, July 2018 experienced much more rainfall (54.8 mm) over 19 events, compared to July 2017. While the frequency of daily rainfall events was greater during August 2018 (n = 17) compared to August 2017 (n = 13), the total cumulative rainfall was slightly less during August 2018 (75.8 mm). The most intense rain events during August 2017 (5.3–16.8 mm) occurred over multiple day periods compared to the most intense rainfall events during August 2018 (6.8–15.4 mm), which were more dispersed throughout the month. Notably, 2017 saw 27.7 mm and 37.6 mm of cumulative rainfall during 13–14 August and 26–30 August (2017), respectively.

Temperature was also variable seasonally and interannually. While peak 7 day mean temperatures consistently reached close to 20 °C during each summer, 2016 and 2018 experienced cooler spring-summer months than 2017. June 2016 and 2018 mean monthly temperatures were 11.5 and 11.9 °C, respectively compared to 13.0 °C for June 2017. July 2016 and 2018 experienced mean month temperatures of 14.9 and 16.0 °C, respectively compared to 17.0 °C during July 2017. August 2018 experience much cooler temperatures (mean monthly = 8.8 °C) compared to August 2016 (11.7 °C) and 2017 (12.4 °C).

2. Materials and Methods

2.1. Remotely Piloted Aircraft System Surveys

High-resolution remotely sensed data were acquired using RPAS at Slump1 six times during the period 2016–2019 to provide insight of geomorphological properties and change. Monitoring of Slump1 began 29-July-2016, approximately one month following the initial evacuation of sediment from the site. Subsequent RPAS surveys were conducted during 6-June and 31-August of 2017 and 7-June and 29-August 2018. One more survey was conducted during 6-June-2019. Three different RPAS vehicles were used beginning with a 3DR Solo equipped with a GoPro 4 camera during the period 2016–2017, followed by an eBee+ (2018–2019) and DJI Phantom 4 equipped with a FC6310 camera (2019). The flight paths were all pre-programmed using Pix4D Capture for the 3DR Solo, Sensefly eMotion (version 3) for the eBee+, and DJI GS Pro for the DJI Phantom 4. Flight paths were oriented perpendicular to dominant wind direction for each survey. The eBee+ was flown with a (SODA) RGB camera, a (Parrot Sequoia) multisensor, and a (Sensfly ThermoMap) thermal sensor. Thermal imagery was only acquired 6-June-2019. The Sequoia multisensor included a 16-megapixel RGB camera, and 1.2 megapixel green, red, red-edge and near infrared (NIR) bands (spanning average wavelengths of 550, 660, 735, and 790 nm, respectively).

Georeferencing of acquired images required use of ground control points and post-processing kinematic (PPK) corrections. Ground control points (GCPs) were used during the initial flights with the Solo 3DR kit. However, these positions were acquired during 2016 with a simple handheld GPS unit with relatively low accuracy. On the other hand, images acquired during subsequent surveys (June 2017 and 2018–2019) utilized highly accurate differential Global Navigation of Satellite System (GNSS) survey control established using Spectra Precision SP80 receivers and PPK corrections, which can achieve horizontal and vertical accuracies of 3 and 5 mm, respectively. A spatially-accurate orthomosaic generated for June 2017 was used to improve the accuracy of the 2016 and August 2017 RPAS products by using common features (e.g., fallen trees) to establish GCPs. Flight path and photo georeferencing during the Ebee+ RPAS surveys (2018–2019) utilized PPK corrections using SenseFly’s eMotion (version 3) flight processing software. Corrections for the aircraft positions were applied following flights using rinex output from a base station GNSS survey (within 1 km of Slump1), which collected satellite locational information in tandem with the flights. A Spectra Precision (SP80) receiver was used for the base station survey as well as the previous ≥7 h static surveys that were used to derive the control point position. Corrections for the GNSS static survey coordinates were provided by the Natural Resource Canada’s Precise Point Positioning service (https://webapp.geod.nrcan.gc.ca/geod/tools-outils/ppp.php). All images from the 3DR and Ebee+ surveys were processed using Pix4D photogrammetry software (versions 3-4.4.12). All acquired images overlapped by at least 65% to provide sufficient spatial reference points for photogrammetric processing. Shortwave radiance values were recorded using the Sequoia system and used in conjunction with image calibration charts to correct multispectral (green, red, red-edge, and NIR) reflectance values in Pix4D. Corrected multispectral bands were used to generate green, red, red-edge, and NIR reflectance layers. The spatial resolution of the orthomosaics and digital surface models (DSMs) were 3.25 cm for the 3DR/Hero 4 RPAS and 2.90 cm for the Sensfly SODA sensor (

Table A1. Details of RPAS photo acquisition and photogrammetry processing results.

Date	Platform	Sensor	Survey Area (ha)	Number of Images	Resolution of Output Layer (s) (cm)	Mean Point Cloud Density (per m3)
29-July-2016	3DR Solo	GoPro 4	6.75	144	3.25	114.47
6-June-2017	3DR Solo	GoPro 4	11.15	278	3.25	197.84
31-August-2017	3DR Solo	GoPro 4	11.04	156	3.67	96.45
7-June-2018	eBee+	SODA	53.99	323	2.90	1448.07
29-August-2018	eBee+	SODA	53.79	342	2.91	940.24
6-June-2019	eBee+	SODA	116.44	491	3.00	334.12
		Sequoia	132.12	488/band	12.77	N/A
		Thermo Map	137.68	3561	18.43	N/A
	DJI Phantom 4	Phantom (FC6310 RGB)	2.45	410	0.84	5572.43

). Reflectance layers generated from the Parrot Sequoia multispectral sensor and SenseFly ThermoMap thermal sensor yielded 12.8 and 18.4 cm pixel resolution, respectively.

The RPAS flight using the Phantom 4 with the FC6310 camera was conducted at a lower altitude (30 m above the top of the slump headwall) to provide higher-resolution mapping data products. The Phantom 4 being used was customized by UAV-Design (Finland) to include GNSS recording capabilities and a spatially-calibrated camera. ReachView flight planning software was used in the field, while Pix4D photogrammetry software was used to generate final mapping products. The Phantom 4 survey also including use of 10 GCPs (five surrounding the slump and five positioned within in the slump) that were corrected using survey control. The spatial resolution of the RGB-derived orthomosaic and DSM was 0.84 cm.

2.2. RPAS Data Analysis

A point cloud of Slump1 and surrounding area was generated using SfM in Pix4D for each sampling interval. The pix4D fill volume tool was used to calculate the cumulative volume of sediment exported from the thaw slump prior to each survey. A polygon extending across the plane of the thaw slump and its entire perimeter provided the elevation reference for determining the volume of space missing from the ground below (i.e., fill volume). An estimation of error was identified for the fill volume of each survey interval using approaches described by Pix4D [34], which assumes that uncertainty is within 1.5 times the grid sample distance (GSD or resolution) volume for each pixel (i.e., 1.5 × GSD³).

A ground irradiance layer was generated, representing the total amount of irradiance (W·m⁻²) the ground surface was exposed to during daylight hours (10:00 to 18:00 local time). The algorithm was implemented using the r.sun.hourly add-on and r.sun modules for GRASS GIS software (version 7.8.3; [35,36]). The required spatial data layers included the RPAS-derived DSM, and slope and aspect layers, which were generated from the DSM. Other input parameters included date and time range. The timing of summer solstice (20-June-2019) was used to model the most intense irradiance values, when the sun height (vertical solar angle) was greatest.

RPAS mapping products were used to explore geomorphological characteristics within Slump1. Differences in vertical profiles and headwall retreat were compared among ten transects that extend from their intersection at the central outflow of Slump1, outward to 5 m past the 2019 headwall. The orientation of the transects were set at 12.5–20° intervals, which correspond with sun azimuth experienced between 6:55 and 17:45 during summer solstice (20-June-2019). The aim was to compare the influence that varying amounts of incident radiation and ground properties (NDVI and temperature) had on thaw slump characteristics and change. For 2019, multispectral (green, red, red-edge and NIR), thermal and irradiance layers were used to classify geomorphological terrain types for Slump1. Maximum likelihood (supervised) classification was used to distinguish terrain types. Pixel training data were identified using the high-resolution (0.84 cm) RPAS-derived orthomosaic and in situ oblique photographs acquired during field campaigns. Accuracy of the supervised classification was evaluated using calculated metrics including overall accuracy, user and producer’s accuracy, kappa statistic, and a confusion matrix output from GRASS GIS [37]. Surface water flow pathways across the varying terrain types were identified using a combination of two terrain modelling procedures and manual digitizing of flow channels visible in the 0.84 cm orthomosaic. The terrain modelling included the topographic water index (TWI), and flow accumulation, both of which were calculated using the DSM for Slump1. TWI was derived using the r.topidx module in GRASS GIS, which utilizes upstream hillslope and area to derive the amount of potential water that drains through each pixel [38,39]. The utility of surface water flow pathway modelling was also explored using the r.flow module for GRASS GIS, which provided a measure of flow accumulation based on density of flowlines passing through each pixel [40].

2.3. Sediment Samples

A total of 40 (thawed) sediment samples were collected from inside Slump1 to characterize physical and chemical properties. During July 2016, six samples were collected along a vertical profile of the northwestern headwall at 10, 60, 80, 100, 120, and 140 cm below the uppermost peat layer horizon. Six more samples were obtained during June 2017 from random locations within low-lying wet areas of the thaw slump floor. During June 2019, samples were collected 10 m apart along three horizontal transects (15 m apart) spanning the floor of Slump1. Sample locations along the transects were recorded using differential GPS. A qualitative assessment of soil moisture (i.e., wet, medium wet, and dry) was recorded at each sampling location along the transects. Approximately 200 cm³ of thawed sediment were collected at each sampling location using cleaned plastic tubing. Once excess sediment was scraped off, the 4 cm tubing sections were capped with minimal air inside, wrapped with plastic film and sealed with tape. Samples were stored at 4 °C until processed at the Water and Environment Lab at Brock University. The samples were analyzed for loss on ignition (LOI), bulk density, gravimetric and volumetric water content, and bulk organic carbon (SOC) and total nitrogen elemental compositions.

Bulk dry soil organic carbon (SOC) and total nitrogen (TN) elemental compositions were measured using methods described by Wolfe et al. [41]. The sediment samples were rinsed with 10% hydrochloric acid to remove carbonates and rinsed repeatedly with de-ionized water until neutral pH returned. Acid-washed samples were freeze-dried to remove moisture and a 500 µm sieve was used to remove coarse organic material. Fine fraction samples (<500 µm) were then analyzed for organic carbon and nitrogen elemental composition using a continuous flow isotope ratio mass spectrometer (CF-IRMS) at the University of Waterloo Environmental Isotope Laboratory.

Relations among sediment analysis results (LOI, dry bulk density, SOC, and TN) and RPAS-derived terrain characteristics were explored to evaluate the influence of thaw slump geomorphological influence on sediment properties. Dry bulk density, SOC, and TN data were also used in conjunction with the cumulative fill (sediment export) volume for each survey interval to estimate SOC and TN mass exported to the debris tongue. The dry bulk density values were converted to mass per unit volume (kg m⁻³) and multiplied by SOC and TN concentrations. These values were multiplied by Slump1 total fill volume to provide estimates of the cumulative SOC and TN mass exported.

3. Results

3.1. Geomorphology of Slump1

Digital elevation models derived from SfM were used to calculate the total volume of sediment exported at the time of each RPAS survey (Figure 3,

Table 1. Slump1 volume and areal increases during the 2016–2019 survey campaigns.

Date	Total Volume (m3)	Accuracy +/− (m3)	Volume Increase (%)	Total Area (m2)	Area Increase (%)
29-July-2016	29,174	295		6342
6-June-2017	32,054	345	10	7163	13
31-August-2017	40,290	477	26	8806	23
7-June-2018	40,637	393	1	9006	2
29-August-2018	47,509	445	17	10,267	14
7-June-2019	48,019	130	1	10,391	1

). By 29-July-2016, approximately one month following the initial evacuation of Slump1, a total 29,174 m³ was exported to the debris tongue, where it was available for transport by the Old Crow River. The exported volume and total area grew by 10% and 13%, respectively, by the following spring (6-June-2017). During the next 86 days (until 31-August-2017), Slump1 experienced the greatest volumetric (26%) and areal (23%) increases during this study. This interval was characterized by the warmest temperatures and average rainfall. Increases in thaw slump volume and area were minimal (1% and 2%, respectively) from 31-August-2017 to the following spring (7-June-2018). During the following 83 days (until 29-August-2018), which were characterized by cooler and wetter conditions than during summer 2017, Slump1 experienced a 17% increase in volume and 14% increase in area. The volume and area of Slump1, again, only experienced slight (1%) increases during the next interval up to the final RPAS survey (6-June-2019). The final cumulative volume of sediment exported by this time was 48,019 m³.

The data acquired using RPAS during 6-June-2019 was used to generate several mapping products (DSM, irradiance, NDVI, and thermal), which were used to identify key geomorphological characteristics of Slump1 (Figure 4). The transects within Figure 4 provide a reference of sun altitude (denoted by colour) and azimuth (line direction) during summer solstice when radiative intensity was at maximum potential. The DSM (hillshade) effectively shows the variable relief within the RTS. This is complemented by the total irradiance (11:30–17:30 UTC-7) layer, which identifies zones within the RTS that received high versus low amounts of incident radiation on clear-sky days. These areas were located at the high-relief locations on the RTS floor, and particularly, towards the northern headwall closer to the ends of transects 6 and 7, where ice-wedge cross sections were often observed in the headwall exposure. These transects were also aligned with the sun at peak altitudes (44.7 and 45.5 °C, respectively) during summer solstice. The more shaded areas were situated immediately northwest of raised areas on the RTS floor, and especially, north of the southwest boundary of the RTS. NDVI shows the lack of vegetation across the majority of the RTS. NDVI values were higher where there was vegetation including ‘older’ black spruce trees (with the highest NDVI values) and moss mat, which had been transported toward the central area of the RTS floor during evacuation events. The relatively high NDVI values within the shaded area in the southwest zone highlight newly established graminoids, which were observed for the first time growing during the 6-June-2019 survey. The thermal data are highly variable. Maximum temperatures (up to 48 °C) are along the steep shoreline of the Old Crow River, facing approximately 150° SSW. Other areas with high temperatures are the raised areas in the eastern and northwestern portions of the RTS floor. These areas had high irradiance, although the other central area with high irradiance had relatively lower surface temperatures. The moisture from the thawing northern headwall kept the sediment cool along these low-lying areas that extended from the northern headwall down to the mouth. The shaded area on the south side also had relatively cool temperatures.

Differences in headwall retreat length during the period 2016–2019 (after the initial evacuation had occurred) were compared among 10 transects extending from the central outflow of Slump1 to 5 m outside of the 2019 headwall boundary (Figure 5;

Table 2. Summary information for the transect profiles 1–10 within Slump1. The timing of solar alignment identifies when the profiles were aligned with the sun’s azimuth during summer solstice (20-June-2019). Mean values along transect profiles are indicated for each RPAS raster layer as well as the total headwall retreat length (m) that occurred during the period 2016–2019. The profiles that experienced the greatest headwall retreat (profiles 5–7) are highlighted in grey.

	Solar Details (20-June-2019)			Mean Values Along Profiles (2019)
Profile	Sun Azimuth Alignment	Sun Altitude	Time	NDVI	Temperature (°C)	Total Irradiance (W·m−2) 11:30–17:30, 20-June	Headwall Retreat 2016–2019 (m)
1	61.5	14	6:55	0.34	18.1	2728	11.7
2	81.5	22	8:25	0.25	17.2	2937	10.4
3	101.5	30	9:50	0.17	17.9	3318	17.6
3.5	111.5	33.6	10:30	0.13	20.5	4155	24
4	121.5	36.8	11:08	0.11	20.6	4123	16.9
5	141.5	41.8	12:20	0.12	19.6	3891	24.9
6	161.5	44.7	13:25	0.13	20.0	4074	27.3
7	181.5	45.5	14:26	0.14	21.4	4810	29
8	201.5	44.4	15:28	0.13	21.9	4643	14.9
9	221.5	41.2	16:32	0.23	21.2	4194	9.9
10	241.5	35.9	17:45	0.27	20.3	4628	1.6

). A transect following the outflow path from the central area (transect 11) is also included in the plot with transect 6. For all ten transects combined, the mean headwall retreat from 29-July-2016 to 6-June-2019 was 17.1 m. The minimum retreat occurred closer to the outflow mouth along SW-NE transects (1, 2, and 10). Greater headwall retreat occurred towards the northwest represented by transects 3, 3.5, and 4, which experienced 17.6, 24.0, and 16.9 m cumulative headwall erosion, respectively. Transects 8–10 are directed to the northeast into an adjacent RTS scar. The greatest retreat during the three-year period was experienced along the northern headwall among transects 5, 6, and 7, which showed 24.9, 27.3 and 29.0 m cumulative headwall erosion length, respectively. Notably, the greatest amount of headwall retreat occurred along transect 7, which aligns with the sun azimuth (181.5°) at the time of maximum sun altitude (45.5° at 14:26 UTC-7; Table 2) during summer solstice (20-June-2019).

RPAS-derived data from 2019 (Figure 4) were evaluated along the profiles (superimposed on RPAS data, e.g., Figure 4a) to compare ground conditions in the directions of varying rates of headwall retreat and sun azimuth and altitudes (Figure 6, Table 2). The majority of Slump1 had no vegetation ground cover as represented by transects 3–8, which have low mean NDVI. The spikes in NDVI along the profiles represent the older vegetation (trees, shrubs and moss mats) that have been transported into the slump during ablation and evacuations. Interestingly, some of this vegetation was relatively healthy during field campaigns (e.g., upright spruce trees) and has potential to re-establish within the RTS scar if evacuations cease. The transects leading to the lateral areas within the southern portion of Slump1 (1, 2, 9 and 10) had the greatest mean NDVI values because of an abundance of older, pre-slump vegetation as well as newly established graminoids. The grasses and mosses present in the stationary sediments among transects 1–3 were observed for the first time during the final survey campaign, almost three years following the initial evacuation. The southwestern area experienced the least amount of headwall retreat and received the greatest amount of shading as shown by the relatively low temperatures and low total irradiance along transects 1 (Figure 4 and Figure 6c,d). There were also relatively low temperatures along the central floor of Slump1 (e.g., along portions of transects 5–7), despite high irradiance (Figure 6d). This area was barren of new vegetation. Field observations indicated that these areas had relatively high moisture supplied from the thawing ice-rich northern headwall. Deep-seated mass wasting was observed along these low-lying corridors immediately following the initial survey (29-July-2016). The sediment with within areas of relatively high ground temperatures (e.g., along transects 4 and 8) had low NDVI, and high irradiance, were dry, and raised (Figure 6a). These raised lobes comprised mostly clay and were less prone to movement.

In summary, the combination of irradiance, multispectral (calibrated green, red, red-edge, and NIR reflectance bands) and thermal data distinguished key geomorphological characteristics within Slump1 (Figure 7), which were observed in the field and were discernable within the high (0.84 cm) resolution orthomosaic. Four classes of interest identified in the field and in the high-resolution RPAS data can be described as (1) low-lying, mixed clay-silt, and barren sediment with high irradiance, low temperature, and prone to mass wasting; (2) headwall-shaded areas with varied relief, cool temperatures and establishment of graminoid and moss vegetation; (3) raised, barren, dry clay lobes, with high irradiance, and temperature; and (4) pre-thaw slump vegetation (e.g., spruce, shrub, grass and moss mat). Figure 7 shows example locations of these terrain classifications identified within Slump1 during the 2019 survey (Figure 7a). The associated normalized spectral signatures of the multispectral, thermal, and irradiance bands within the polygons (Figure 7b) have super-imposed circles identifying the key differences that allow for effective differentiation of the terrain types. The visible and NIR bands were useful for distinguishing all four classes, but mostly for identifying pre-slump vegetation and shaded areas with new vegetation growth (Figure 7b,c). High irradiance values distinguished the open and barren terrain types (classes 1 and 3) from the shaded areas (class 2; yellow circle in Figure 7b,c). Whereas thermal data effectively distinguished the warmer, raised dry clay lobes (class 3) from the cooler, low-lying high flow class 1.

To illustrate further, these RPAS-derived mapping layers from the 2019 survey were used to generate a terrain classification map to highlight the geomorphological differences within Slump1, (Figure 8). The maximum likelihood classification approach was implemented using the training data identified in Figure 7a in GRASS GIS software (version 7.8.3). The overall accuracy of the classification was 98.2%, with a kappa coefficient of 0.97. The most confusion among classes was between shaded, moist areas and the open, moist areas (user accuracy = 90.2%). The other groupings had user and producer accuracies between 97.4 and 99.9%. Terrain classification was attempted for previous years using only RGB and DSM (irradiance) layers with mixed success. While the general areas of the terrain classifications spatially aligned with the 2019 products, user and producer accuracies decreased to as low as 50%. Therefore, the thermal and NIR reflectance data acquired during 2019 provided key layers for differentiating thaw slump geomorphic characteristics.

Surface water flow pathways within Slump1 were identified using the combination of modelling approaches (TWI and flow accumulation) using the DSM and manual digitizing of visible flow channels identified in the 0.84 cm resolution orthomosaic. Based on qualitative comparison of the modelled results with the orthomosaic, TWI (values > 4) provided favorable identification of relatively lower-flow pathways. However, the TWI did not detect much of the larger, incised channels, mostly located in the low-lying areas of the thaw slump. The sensitivity of TWI in these areas was likely reduced because of the effects of shadows within the RPAS-derived images that were used to generate the DSM. Therefore, manually digitized channels (n = 196) were required to enhance identification of the high surface water flow pathways. The flow accumulation model (from the r.flow module in GRASS GIS, version 7.8.3) provided favorable results for identifying the raised ridges. For our application, these mapping layers highlighted differences in low- and high-flow water accumulation among classified thaw slump terrain types (Figure 8). Locations of ridges along the raised dry clay lobes are highlighted in red, which are typically perpendicular to the headwall. Surface water flow pathways (highlighted in blue) represent areas that convey overland water flow and the suspended sediment eroded from fluvial processes. Therefore, many of these high liquid-flow channels are incised into the thaw slump floor.

Degradation and aggradation differences among terrain classification types were compared for each time interval between survey campaigns (Figure 9). The pre-slump vegetation class was not included in the analysis. Headwall degradation (<−2 m), which provided the greatest source of sediment for transport along thaw slump flow paths was also not included since the focus was to compare differences among terrain classes defined within the area of the thaw slump floor. Pixel sampling zones of each class for each survey were identified using the orthomosaic and in situ photographs acquired during the beginning date of each interval. Results were plotted as scaled density for semi-quantitative assessment rather than frequency since variable spatial portions of the terrain type areas were used (Figure 9). Degradation was dominant during the first two intervals up to August 2017 for the three terrain classes (~−0.25–−1.25 m). During this interval, the shaded southwestern zone also experienced some aggradation in spots (<1 m). Aggradation (between 0.2 and 0.9 m) was experienced during fall 2017 to spring 2018, mostly for the low-lying, high-flow accumulation pathways. This sediment was exported during the period June–August 2018. During this interval, similar amounts of erosion occurred within the raised dry clay zones, while less degradation was experienced within the shaded southwest zone. Elevation change was negligible for the three terrain classification types during fall 2018 to spring 2019. Overall, the low-lying, high-flow accumulation pathways experienced the most dynamic elevation changes since these zones received sediment from the headwall as well as the raised dry clay zones. The raised dry clay zones generally degrade over time since there is no sediment source for aggradation. The thawing headwall sediment accumulates within the lower-lying areas that are more prone to evacuation during mass wasting events. The shaded southwest zone was more stabilized than the other terrain types during the last year of this study.

3.2. Thaw Slump Sediment Organic Carbon, and Total Nitrogen

Percent SOC and TN values show variability when grouped by terrain unit classifications (Figure 10,

Table A2. Sediment chemistry and RPAS data values for samples and sample locations. RPAS multispectral (green, red, re-edge, and near infrared) reflectance data and DSM-derivative (slope and aspect) data are represented by majority (modal) pixel values within 1 m radius zones around the transect points (n = 28) sampled during 6-June-2019. Terrain ID values represent low-lying high flow (1), shaded with new vegetation (2), and raised open dry clay lobes (3).

Sample ID	Sample Date	Dry Bulk Density (g/cm3)	SOC (%)	TN (%)	Terrain ID	Slope	Aspect	Irradiance (W m−2)	Green	Red	Red-Edge	NIR	Thermal (°C)
SED1	6 June 2019	0.871	1.980	0.173	2	16.383	45.0	1365.804	0.0610	0.0699	0.0786	0.0903	15.785
SED2	6 June 2019	1.189	1.830	0.128	2	0.625	45.0	2065.371	0.0580	0.0652	0.0778	0.0943	15.644
SED3	6 June 2019	0.558	0.910	0.085	3	1.171	315.0	4252.25	0.0821	0.0899	0.0951	0.1114	19.839
SED4	6 June 2019	0.728	1.100	0.090	3	11.735	360.0	3550.799	0.1046	0.1204	0.1203	0.1382	17.745
SED5	6 June 2019	0.795	0.990	0.091	1	2.890	360.0	3919.699	0.0687	0.0752	0.0801	0.0981	16.955
SED6	6 June 2019	1.240	1.430	0.133	1	0.000	270.0	4246.106	0.0591	0.0595	0.0675	0.0869	18.514
SED7	6 June 2019	1.076	1.310	0.114	1	7.479	315.0	4414.696	0.0612	0.0693	0.0857	0.1069	17.918
SED8	6 June 2019	1.089	0.790	0.072	1	0.103	270.0	4400.13	0.0550	0.0613	0.0767	0.0911	18.856
SED9	6 June 2019	1.100	0.770	0.070	1	5.030	251.6	4688.197	0.0826	0.0986	0.1219	0.1352	20.591
SED10	6 June 2019	1.171	0.890	0.081	3	5.622	135.0	1683.798	0.0486	0.0564	0.0640	0.0766	13.393
SED11	6 June 2019	1.063	3.240	0.215	2	20.158	90.0	2945.124	0.0511	0.0568	0.0723	0.0945	13.465
SED12	6 June 2019	1.249	1.780	0.153	2	15.516	45.0	1945.984	0.0419	0.0387	0.0532	0.0732	14.755
SED13	6 June 2019	1.235	1.550	0.111	1	3.656	315.0	3700.076	0.0584	0.0641	0.0739	0.0879	16.513
SED14	6 June 2019	1.142	1.510	0.137	3	10.666	296.7	3412.757	0.0893	0.1026	0.1132	0.1268	19.397
SED15	6 June 2019	1.222	1.580	0.148	1	3.276	315.0	4180.191	0.0776	0.0854	0.0927	0.1081	17.732
SED16	6 June 2019	1.312	2.450	0.169	1	0.000	315.0	4245.879	0.0689	0.0753	0.0950	0.1077	17.491
SED17	6 June 2019	1.156	1.000	0.085	3	3.280	315.0	3327.682	0.0607	0.0686	0.0763	0.0871	19.482
SED18	6 June 2019	1.074	1.600	0.136	3	1.045	225.0	5019.789	0.0914	0.1040	0.1149	0.1286	20.005
SED19	6 June 2019	1.177	1.430	0.115	3	7.242	315.0	4416.411	0.0850	0.1113	0.1178	0.1327	21.317
SED20	6 June 2019	1.200	1.240	0.096	1	1.397	225.0	3666.532	0.0566	0.0689	0.0803	0.0957	18.215
SED21	6 June 2019	0.968	1.900	0.127	1	2.093	270.0	4233.487	0.0495	0.0562	0.0763	0.0977	22.304
SED22	6 June 2019	1.183	2.960	0.181	2	12.832	45.0	405.6073	0.0365	0.0321	0.0474	0.0636	14.379
SED23	6 June 2019	1.283	0.960	0.086	1	0.640	315.0	3334.005	0.0666	0.0695	0.0836	0.0980	15.444
SED24	6 June 2019	0.808	1.700	0.127	3	14.883	315.0	4886.898	0.0799	0.0873	0.0988	0.1198	19.672
SED25	6 June 2019	1.143	1.860	0.182	1	5.843	270.0	4527.966	0.0633	0.0682	0.0764	0.0890	15.073
SED26	6 June 2019	1.264	0.470	0.045	1	0.000	0.0	4245.756	0.0692	0.0750	0.0898	0.1047	12.915
SED27	6 June 2019	1.286	1.560	0.124	1	0.478	225.0	4271.771	0.0309	0.0222	0.0320	0.0428	18.681
SED28	6 June 2019	1.158	1.810	0.135	3	12.537	270.0	4204.556	0.1023	0.1244	0.1429	0.1571	20.549
exposure_10cm	29 July 2016	1.606	0.955	0.063	1
exposure_60cm	29 July 2016	1.585	0.781	0.070	1
exposure_80cm	29 July 2016	1.457	0.624	0.060	1
exposure_100cm	29 July 2016	1.471	0.761	0.070	1
exposure_120cm	29 July 2016	1.604	0.850	0.073	1
exposure_140cm	29 July 2016	1.506	0.794	0.065	1
west1	7 June 2017	1.474	0.992	0.081	1
south1	7 June 2017	1.597	0.921	0.076	1
south2	7 June 2017	1.462	1.587	0.122	1
west2	7 June 2017	1.377	1.139	0.078	1
south3	7 June 2017	1.474	1.106	0.093	1
unsure	7 June 2017	1.299	2.496	0.161	1

). Although the number of samples obtained from each terrain grouping are unequal, there are notable distinctions. The greatest differences in SOC and TN concentrations were identified between the shaded areas (with grasses and moss) with the highest values and the raised dry clay lobes with lower values. Twenty-two of the 26 samples from the low-lying, high-flow accumulation zone had SOC values less than the minimum value (1.78%) in the shaded southwest zone. The majority (8/9) of samples from the raised, dry clay lobes also had SOC less than the minimum value from the shaded area. Similarly, TN (%) was typically lower for samples outside of the shaded southwest zone. Thawed sediment sampled from the northwester headwall during July 2016 at 10–140 cm below the ground surface also had relatively low SOC and TN values. Twenty-one of the 26 samples from low-lying, high flow zones and six of the nine samples from raised dry clay zones had TN (%) less than the minimum (0.13%) value in the southwest shaded zone. The low-lying, high-flow accumulation areas had the most variable SOC and TN concentrations since they convey sediment from the other geomorphic zones (headwall, shaded southwest area, and raised dry clay lobes).

More detailed exploration of the relations among terrain and sediment characteristics is shown in Figure A1. Principal components analysis (PCA; Figure A1a) revealed that the higher SOC and TN concentrations in the shaded terrain type was associated with greater slope, NDVI, and bulk density, and lower aspect, temperature, and irradiance. On the other hand, sediment from the raised dry clay lobes was strongly associated with higher aspect, temperature, and irradiance. Samples from the low-lying high flow zones plot between the other terrain classes. Scatterplots of SOC and TN with the other terrain characteristics (Figure A1b) highlight the distinction of the shaded terrain class from the other classification groupings based on individual RPAS-derived datasets. This is most notable for slope, aspect, irradiance and temperature. TWI was not distinguishable among terrain types. There is much more overlap between the low-lying high flow and raised dry clay lobe classifications, especially in terms of nutrient concentrations. However, there is notable distinction in slope and temperature values between these two groupings with the raised dry clay lobes typically being steeper and warmer.

Sediment analysis was coupled with RPAS-derived volume calculations to estimate the total amount of SOC and TN exported as sediment from Slump1 between each sampling interval. Of the 40 sediment samples (Table A2), the 26 samples that were collected within low-lying, high sediment flow (evacuation-prone) areas of the thaw slump were used to calculate mass of SOC and TN export to the debris tongue on the Old Crow River. These samples included 14 from along the 2019 transects and 12 from non-transect locations collected prior to 2019. Median dry bulk density, and median SOC and TN concentrations were 1.293 g cm⁻³, 1.05% and 0.09%, respectively (Figure 10). SOC and TN storage estimates in the high-flow areas of Slump1 were derived from multiplying bulk density (converted to kg m⁻³) by SOC and TN concentrations. Median SOC and TN values were 14.9 kg C m⁻³, and 1.2 kg TN m⁻³, respectively. After multiplying these values by fill volume values from each RPAS survey, we estimated that Slump1 exported approximately 713 t SOC and 58 t TN to the Old Crow River from the time of initiation during June 2016 to June 2019. Following the initial detachment failure during June 2016, the most SOC (122 t) and TN (10 t) was exported during summer 2017 compared to during 2018 (SOC = 102 t and TN = 8 t).

4. Discussion

This study demonstrates a unique approach that integrates RPAS-derived datasets and in situ measurements for monitoring thaw slump geomorphology, drivers, and influence on soil properties and nutrient mobilization. Discussion of this multifaceted study is broken down among volumetric change monitoring, associated drivers (meteorology, orientation), thaw slump geomorphology, and detected relations with soil nutrients, which can be utilized to refine estimates of downstream impacts. Our findings underscore the dynamic geomorphic properties of permafrost retrogressive thaw slumps that should be considered when anticipating associated impacts on important Arctic and subarctic landscapes such as OCF, that are sensitive to ongoing climate change.

The 2016–2019 time series of high-resolution DSMs produced from SfM provide insight of the seasonal and interannual timing of evacuations. This approach has been performed for many sites across Northwest Territories, Canada, where thaw slumps have reached fill volumes of 3.2 × 10³ to 5.9 × 10⁶ m³ [17]. Slump1 reached a fill volume of 4.8 × 10⁴ m³ by spring 2019, which was 61% greater than its size following the first RPAS survey in July 2016. The (~1-ha) area of Slump1 is relatively large for OCF based on comparison with other thaw slump scars observed in oblique aerial photographs from reconnaissance flights and general assessment of airborne datasets (i.e., AVIRIS-NG). Thaw slump area is much more variable in other landscapes. For example, the median area of thaw slumps across Tuktoyaktuk, NT is 0.4 ha, but reach as large as 5 ha [13].

Volumetric change coupled with analysis of meteorological data provide insight of the drivers of Slump1. While there was greater total accumulated rainfall during 2018, greater evacuation of sediment occurred during the warmer 2017 summer when precipitation was more focused within shorter (2–3 day) time periods during August. These findings are consistent with other studies in northwestern Canada [3,7,13], Alaska [42], and Russia [13], where debris often accumulated on thaw slump floors following headwall erosion and was subsequently exported during precipitation events later in the summer. The high-resolution DSMs produced from RPAS-derived images were ideal for evaluating the timing of steady versus mass-wasting events. However, future studies should attempt more frequent acquisition of high-resolution thaw slump size data to refine insight of meteoric influence on processes driving volumetric change. For example, the timing of RPAS surveys over multiple thaw slumps could be aligned with the frequent (11-day interval) acquisition of TanDEM-X satellite data, which was used by Zwieback et al. [13] for evaluating the timing of thaw slump volumetric change across vast landscapes in northwestern Canada and Russia.

The orientation of thaw slumps can be variable according to local geomorphic properties. While local geomorphic drivers including meteorological conditions and erosion along rivers and coastal areas often represent dominant controls on initial detachment failures, the direction of subsequent maximum headwall retreat and ablation can often be linked with orientation. As noted by Lewkowicz [7], the initiation of thaw slumps on Banks Island were triggered by coastal erosion with high incident radiation on south- and southwest-facing slopes, promoting headwall retreat. Pronounced northward headwall retreat has been observed at other notable sites including the Batagay mega slump in northern Yakutia, Russia [16] and in northwestern Canada [10,17]. Our results are consistent with these observations and highlight the utility of modelled (daily cumulative) irradiance (i.e., during 20-June; summer solstice) for highlighting the importance of orientation in driving greater northward headwall retreat. Other studies have indicated that orientation may be less influential on rates of headwall retreat. For example, Swanson and Nolan [18] noted that the growth of thaw slumps in the Noatak Valley, Alaska was not strongly linked with orientation, particularly with high-irradiance angles, with the greatest headwall retreat rates located along west- and north-facing slopes. This can also be the case for coastal environments, where thaw slump activity is perpetuated strongly by coastal thermo-erosion [15,43]. However, there are many landscapes where the orientation of thaw slump activity may be more effectively predicted when incorporating high-resolution DSM-derived irradiance data. While the orientation of hillslope thaw slumps in the Richardson Mountains-Peel Plateau was predominantly east facing [11] and in the Mackenzie Delta predominantly south and west facing [44], the direction of maximum headwall retreat during years following initiation was not identified. As well, Lewkowicz and Way [14] identified drastic increases in the frequency of thaw slumps on Banks Island, with the highest frequency occurring within the southern portion of the study area. Temperature was identified as an important determinant of this general spatial pattern, although it will be interesting to evaluate the influence of orientation on thaw slump growth rates where remotely sensed data of sufficiently high resolution are available.

Acquisition of high-resolution datasets of thaw slumps at varying stages of development will be useful for predictive spatial modelling of future thaw slump activity and hazard mapping. Preliminary assessment of thaw slump spatial distributions in OCF indicates that orientation is an important geomorphic factor of their formation and development. Thaw slumps predominantly exist here on south- and southeast-facing slopes that experience headward erosion along the steep banks of high-order river channels (i.e., Old Crow River and lower Johnson Creek). Future analyses using high-resolution RPAS-derived data and data acquired during NASA’s Arctic-Boreal Monitoring Experiment Airborne Campaigns (https://above.nasa.gov/airborne_2017.html) will refine our inventory of thaw slump frequency and development patterns.

Multispectral, thermal, and DSM-derived datasets provided refined insight of the complex geomorphology of thaw slumps, which influence their long-term development and influence on biogeochemical cycling. Several studies have utilized SfM with RPAS-derived data to monitor thaw slump geomorphology and volumetric changes in northwestern Canada [17,45,46]. Here, we explored the utility of high-resolution multispectral and thermal datasets for classification of thaw slump terrain types that had varying influence on sediment export and biogeochemical cycling. The most dynamic portion of the Slump1 floor was the low-lying, high-flow accumulation zones, which were readily prone to evacuations during rain events. These low-lying areas received sediment from the thawing headwall and raised dry clay lobes. The raised dry clay lobes were less prone to mobility during mass wasting events but experienced more long-term degradation and weathering as sediment eroded from these areas to the lower-lying high-flow accumulation zone during rain events. The raised dry clay lobes comprised clay sediment that were thawed longer than most of the sediment being conveyed seasonally along the low-lying high-flow accumulation zone. The more shaded zone in the southwest area experienced stabilized conditions relatively early, and as a result, saw recent growth of moss and graminoid vegetation. The low-lying, high sediment flow terrain zones were differentiated from the raised dry clay lobes based on temperature and slope. The raised dry clay lobes had high temperature since they were out of the flow pathway of headwall meltwater and, with greater slope, also have more efficient drainage than the low-lying high flow zones. The shaded zone with new vegetation, which also had cooler temperatures, could be distinguished with its unique combination of (green to NIR) spectral signatures and low cumulative daily irradiance.

Past studies have identified the role of fluvial erosion transporting sediment from thaw slumps e.g., [3,38]. Incised channels, formed from fluvial erosion, were observed during field campaigns in Slump1, which at times, directed sediment-laden water to the Old Crow River. Depending on the intensity and duration of rainfall events and ground ice meltwater volumes, effective drainage of runoff along these conduits reduced pooling and supersaturation that may have otherwise promoted deep seated mass wasting. Surface water channels were easily identifiable within the high-resolution orthomosaic. Other flow pathways that were more difficult to visualize and document were effectively modelled using the TWI derived from the DSM. For example, sheetflood pathways, that form alluvial fans at the base of the headwall [8], were effectively identified with TWI (values > 4) and connected with visibly incised channels. Flow accumulation analysis [40] also identified low or non-flow ridges that were less prone to being exported during mass wasting. While shadows can hamper the detection of higher-order incised channels, these high-resolution optical multispectral datasets can provide key insight of the fluvial processes influencing thaw slump development and influence on downstream environments.

At landscape and regional scales, the sediments exposed and mobilized by individual thaw slumps are often characterized by a homogeneous nutrient (i.e., SOC) supply, which varies according to broad regional soil types e.g., [30]. Obu et al. [31] did find that differences in SOC and TN concentrations were associated with feature types including mass-wasting disturbance (i.e., thaw slumps), peat, and fluvial deposits, as well as terrain characteristics including slope angle and TWI. For instance, SOC and TN was found to be positively correlated with TWI, negatively correlated with slope angle and of lower concentrations in the upper sediments of thaw slumps relative to other landscape features. Findings presented here suggest that the influence of slope angle and TWI on SOC and TN can be overshadowed at finer scales by physiographic properties within individual thaw slumps. We detected differences in SOC and TN concentrations among the terrain classifications in Slump1 with the shaded areas having higher values, despite having steeper slope angles. The retention of soil moisture within the shaded area of Slump1 provided ideal microenvironment for establishment of mosses and graminoids, which allowed for accumulation of organic matter. This was contrary to the low nutrient values found in the non-shaded areas that were raised and aligned with the highest incident radiation. These dry clay lobes slowly degraded during this study but retained sediment that had been exposed to the most weathering. Therefore, low SOC and TN are likely in response to leaching downslope, which is consistent with Obu et al. [31]. SOC and TN values from the low-lying high-flow accumulation zones were the most variable since sediment from the other terrain classification zones were mixed and transported along these pathways. Because mixing of sediments occurs during evacuation events, we can expect that outflow sediment nutrient concentrations to be proportional to the size and relative contribution of the source geomorphic terrain types.

The SOC and TN results (14.9 kg C m⁻³ and 1.2 kg TN m⁻³, respectively) used in our nutrient flux calculations are less than those found by Obu et al. [31] on Herschel Island, Yukon, Canada who reported thaw slumps to have 20.9 kg C m⁻³, and 3.7 kg TN m⁻³, respectively, based on core samples. They suggest that the nutrient concentrations were higher at their sites, below 50 cm from the surface, because of greater aeration and microbial activity at the surface [47] and subsequent carbon degradation [48]. This may also be the case for our surface samples from the raised dry clay lobes in Slump1. However, the six samples obtained along a vertical profile (from 10 to 140 cm below surface) at the headwall of Slump1 had similar dry bulk density, SOC, and TN values as the other zones of the thaw slump floor. The consistency of these values suggests that our SOC and TN flux measurements are representative of sediment being exported during mass wasting events. These values are also consistent with results from the North Slope of Alaska [49]. However, additional sediment sample analysis of thaw slumps in OCF should be completed to ensure our calculations are not conservative.

High-resolution terrain classification mapping of thaw slumps using thermal, multispectral and elevation datasets provides useful context when interpreting thaw slump nutrient flux and biogeochemical cycling. Consideration of within-region variation in thermokarst intensity and landscape properties is critical for determining biogeochemical responses [23,25]. Ongoing research should continue to maximize the use of high-resolution elevation, multispectral, and thermal datasets within situ measurements of sediment and downstream properties to build our knowledge of thaw slump biogeochemical cycling and associated impacts on aquatic environments. Widespread high-resolution mapping and analyses can provide much needed context to strengthen estimates of carbon stocks [30], carbon atmospheric flux [27], and carbon mobility projections based on the varying mitigation strategies [50].

5. Conclusions

This study showcases the utility of RPAS datasets for monitoring the activity and geomorphological characteristics of permafrost retrogressive thaw slumps. Furthermore, findings comprise key insight into the processes driving thaw slump development and relations among geomorphic and sediment nutrient compositions. Integrated datasets were used to estimate sediment load and SOC and TN load volumes exported to the Old Crow River between each survey interval during the three-year study.

Structure-from-motion analysis of RPAS-derived RGB photographs provided key information of retrogressive thaw slump geomorphic properties and change in response to meteorological drivers. RPAS surveys revealed that 29,174 m³ of sediment was exported during the initial detachment failure with subsequent evacuations exporting an additional 18,845 m³ until June 2019. Greater sediment export occurred during the warmer and drier summer 2017 (8236 m³) compared to summer 2018 (6872 m³). However, the fewer rain events during August 2017 were of higher intensity and concentrated within shorter time intervals compared to the rain events more evenly distributed during 2018. The northern headwall experienced the greatest amount of retreat owing to its alignment with the highest incident radiation and the ablation of ice wedges that were commonly observed along its face (Figure 4). More intense incident radiation during the warmer summer 2017 accounts for greater headwall retreat and subsequent evacuations, which was enhanced during August rain events. The north-south alignment of the maximum direction of thaw slump headwall retreat showcases the importance of incident radiation on the detected volumetric change. Modelled (summer solstice) cumulative daily solar irradiance (generated from a DSM) is highly useful for explaining anticipating headwall retreat properties.

Combined RPAS-derived irradiance, multispectral and thermal data provided key insight of thaw slump geomorphology. Maximum likelihood supervised classification using these layers provided high-resolution mapping of thaw slump terrain types including (1) low-lying, high-flow accumulation, evacuation-prone sediment; (2) shaded areas with graminoid and moss vegetation; (3) raised, dry clay lobes; and (4) pre-thaw slump vegetation (e.g., spruce, shrub, grass, and moss mat).

Surface water runoff pathways within a thaw slump can be detected using combined modelling and manual visual identification approaches. Small, low-order channels (formed from pluvial and rill erosion) and sheetflood pathways can be readily identified with modelling (e.g., TWI). However, this approach is less effective for highlighting larger channels covered by shadows during the RPAS surveys. On the other hand, these higher-order channels are easily distinguishable with high-resolution (0.84 cm) RGB imagery and can be manually digitized.

SOC and TN concentrations varied among classified terrain types with the highest concentrations within the shaded areas with steeper slope values where new vegetation had established, and the lowest concentrations in the raised dry clay lobes that were subjected to prolonged weathering. The low-lying, high-accumulation zone had the most variable concentrations since these areas conveyed sediment from the other two classes. Coupled RPAS-derived volume, sediment dry bulk density and SOC and TN data indicated that Slump1 exported approximately 713 t SOC and 58 t TN during the three-year study. The most SOC (122 t) and TN (10 t) was exported during summer 2017 compared to during 2018 (C_org = 102 t and TN = 8 t).