Monitoring Inter- and Intra-Seasonal Dynamics of Rapidly Degrading Ice-Rich Permafrost Riverbanks in the Lena Delta with TerraSAR-X Time Series

Abstract

Arctic warming is leading to substantial changes to permafrost including rapid degradation of ice and ice-rich coasts and riverbanks. In this study, we present and evaluate a high spatiotemporal resolution three-year time series of X-Band microwave satellite data from the TerraSAR-X (TSX) satellite to quantify cliff-top erosion (CTE) of an ice-rich permafrost riverbank in the central Lena Delta. We apply a threshold on TSX backscatter images and automatically extract cliff-top lines to derive intra- and inter-annual CTE. In order to examine the drivers of erosion we statistically compare CTE with climatic baseline data using linear mixed models and analysis of variance (ANOVA). Our evaluation of TSX-derived CTE against annual optical-derived CTE and seasonal in situ measurements showed good agreement between all three datasets. We observed continuous erosion from June to September in 2014 and 2015 with no significant seasonality across the thawing season. We found the highest net annual cliff-top erosion of 6.9 m in 2014, in accordance with above-average mean temperatures and thawing degree days as well as low precipitation. We found high net annual erosion and erosion variability in 2015 associated with moderate mean temperatures but above average precipitation. According to linear mixed models, climate parameters alone could not explain intra-seasonal erosional patterns and additional factors such as ground ice content likely drive the observed erosion. Finally, mean backscatter intensity on the cliff surface decreased from −5.29 to −6.69 dB from 2013 to 2015, respectively, likely resulting from changes in surface geometry and properties that could be connected to partial slope stabilization. Overall, we conclude that X-Band backscatter time series can successfully be used to complement optical remote sensing and in situ monitoring of rapid tundra permafrost erosion at riverbanks and coasts by reliably providing information about intra-seasonal dynamics.

1. Introduction

The Arctic is warming almost twice as fast as the global mean [1]. This phenomenon contributes to the rapid degradation of permafrost landscapes which has far-reaching implications for ecological and social systems at local and global scales [2,3,4,5]. Thermokarst and thermal erosion are prominent processes associated with permafrost degradation. Thermokarst is defined as a process by which characteristic landforms result from the thawing of ice-rich permafrost or the melting of massive ice [6], such as thaw slumps [7] and thermokarst lakes or basins [8]. Thermal erosion is the removal of ice-bearing permafrost by the combined thermal and mechanical action of moving water [6,9], and is also the most effective process of permafrost degradation [10,11]. It dominates the reshaping of permafrost landscapes by the erosion of ice-bearing riverbanks [12,13], coastlines [9,14] and lake shores [15] and leads to the formation of retrogressive thaw slumps [7], as well as gullies and valleys [16]. Thermo-denudation is the combined influence of solar insulation and heat advection on the energy balance of the ground surface causing erosion above the water level [9].

Riverbank erosion in continuous permafrost is strongly influenced by the hydro-mechanical and thermal forces of water [9] but the magnitude and timing is determined by factors such as sediment type, snow and ice cover, air and water temperature, and wind [12]. Thermo-erosional niching is the undercutting of cohesive, ice-rich banks by water and is a prominent process of riverbank erosion [17] leading to rates of up to 19 m per year [11]. Different stages of erosion are observed over long time periods with riverbanks eventually stabilizing, a process that can occur in as little as three years [11]. In later stages of riverbank erosion, thermo-denudation in the upper part of the bank is the main driving process leading to bank erosion. Seasonal stages define the relative contributions of thermal and mechanical components of bank retreat [18]. In the lower Lena River, Costard et al. [13] identified a strong seasonality in the erosion of ice-rich riverbanks with flooding and ice push during and after ice break-up in the spring as the main drivers of erosion. High variability in erosion data has been previously observed at eroding ice-rich riverbanks [19] and is possibly linked to varying ground ice contents along bank sequences.

Ice-rich permafrost regions in Siberia have been studied extensively with regard to thermal erosion using optical remote sensing data [20,21]. Yet, the intra-seasonal dynamics of rapid permafrost degradation (such as the erosion of permafrost coasts or riverbanks) remain largely unknown due to a lack of intra-seasonal, high temporal resolution measurements. This is mainly due to persistent cloud coverage in the Arctic and the associated limitation in successful optical image acquisitions. High spatial resolution active microwave remote sensing is a promising alternative as it is unaffected by cloud cover and solar illumination. Synthetic Aperture Radar (SAR) antennas transmit pulses in the microwave length of the electromagnetic spectrum and receive their echoes after they propagate through the atmosphere and are scattered by objects on the earth’s surface. The received microwave signal results from the physical structure and the dielectric properties of the surface in addition to parameters of the satellite viewing geometry. The resulting backscatter intensity can be used to distinguish between different land surfaces.

Today’s SAR systems can provide high spatial resolution data and operate using a variety of wavelengths and polarization configurations depending on the desired application. However, in contrast to optical remote sensing, the interpretation and evaluation of SAR data from highly dynamic tundra environments is a newly developing and challenging area of research. Through joint international programs and efforts by the Polar Space Task Group (PSTG), SAR time series at high spatial and temporal resolution are increasingly available in selected tundra environments [22] offering a unique opportunity for detailed assessment of intra-seasonal permafrost degradation. SAR data has been successfully used to characterize Arctic tundra environments [23,24] and to monitor surface changes [25,26,27]. X-Band SAR data have also been used with differential interferometric SAR techniques to detect and monitor mass movements of landslides and rock falls in a variety of locations and conditions [28,29]. Single pass interferometry with X-band SAR data has been tested for the detection of permafrost degradation features in the Arctic [30,31] and reported challenges in applying interferometric SAR to complex terrain in Arctic environments. Zwieback et al. [31] tested interferometric methods at the same study site presented in this paper and found that the complex topography with foreshortening and layover effects introduced a strong azimuth ambiguity. In addition, Antonova et al. [27] reported low coherence values over tundra surface on Kurungnakh Island, Siberia (where our site is located) using 11 day SAR multi-pass interferometry. Investigations from Widhalm et al. [24], Antonova et al. [27] and Ullmann et al. [32] on X-Band SAR data in tundra landscapes showed low spatial variability in X-Band backscatter for different tundra types (polygonal wet tundra and shrub tundra) during the snow free summer months. In contrast, bare ground conditions that occur after erosion and are accompanied by the removal of the vegetation layer have been shown to have higher backscatter in X-Band with HH polarization [32]. The stability of the tundra backscatter and the high backscatter of bare ground allows for differentiation of these two surfaces.

In this paper, we use a three-year X-Band time series acquired by the TerraSAR-X (TSX) satellite to examine the capability of X-Band SAR to monitor erosion of an ice-rich riverbank in the central Lena Delta, Siberia at high temporal and spatial resolution. Our specific objectives are to (1) evaluate the performance of X-Band SAR data acquired by TSX for the analysis of rapid permafrost erosion against in situ and optical datasets; (2) quantify inter- and intra-annual cliff-top erosion of a rapidly eroding ice-rich permafrost riverbank; (3) investigate relationships between cliff-top erosion and climate data; and (4) examine backscatter intensity dynamics of tundra and cliff surfaces.

2. Study Area

The study site is located within the central Lena River Delta in northeastern Siberia between 72°–74° North and 123°–130° East (Figure 1a). The Lena River Delta is the largest of the Arctic deltas covering about 29,000 km² [33,34] and is located in the zone of continuous permafrost with a maximum permafrost thickness of 500 to 600 m [35]. The region is characterized by an Arctic continental climate with strong seasonal variations. The mean annual air temperature is −12.5 °C, with mean monthly temperatures in February and July of −33.1 °C and 10.1 °C, respectively [36]. Rainfall occurs mostly from mid-May to the end of September with a seasonal mean of 125 mm [36]. Snow cover in the Lena Delta develops between October and July with an average thickness of about 40 cm [36].

We collected erosion data on a 1930 m long rapidly eroding riverbank section on the southeastern side of Kurungnakh Island located in the Olenekskaya channel, which is one of the four main branches of the Lena River Delta (Figure 1b). The sediments at the study site consist of fine-grained, ice- and organic-rich material that is widely distributed in Siberia and is also referred to as Yedoma [37,38]. The distribution of Yedoma sediments in the Lena Delta is shown in Figure 1a. Maximum heights of the undisturbed tundra upland on top of the Yedoma at the study site are about 55 m a.s.l. [39] (Figure 1c). The Yedoma at the study site features polygonal networks of large ice wedges up to 20 m deep and 5 to 7 m wide that can be exposed along the riverbank [40]. The Yedoma is underlain by a 15 to 20-m thick fine layered sand unit of fluvial origin [41]. Gravimetric ice contents of the Yedoma in this area can reach up to 150 weight %, related to the dry sediment weight [42]. The study site is located in a riverine environment without marine wave activity, about 10 km inland from the outer delta margins. The annual average water discharge of the Olenekskaya channel shows strong annual variations ranging from 700 to 1.500 m³/s [43] and open water conditions exist from June to October. The water levels during spring river ice break-up in late May to early June can reach up to 12 m above normal at this location but does not exceed the upper limit of the sand unit (Figure 1c). While the lower sand unit experiences fluvial erosion, the Yedoma cliff top is not directly affected by thermal erosion through Lena River waters, even during spring floods. Thermo-denudation processes at the cliff top are thus mainly driven by atmospheric influences, similar to other coastal Yedoma sites in the Laptev Sea [14], resulting in the formation of thermo-cirques, also referred to as thaw slumps. Thermal erosion is occurring at various cliff locations with complex geometry where gullies channelize runoff from the tundra of the Yedoma hinterland.

The studied cliff section is covered by small mounds of partially thawed and revegetated silty material called Baydzarakhs which remain after the thaw of ice wedge networks (Figure 1c). Steep cliff headwalls feature large ice wedge exposures and indicate very high ground ice contents. Active cliff sections that undergo erosion have bare and wet ground and meltwater streams can result in channel formation. Non-active cliff sections are drier and feature pioneering vegetation. These areas show similarities to the floors of retrogressive thaw slumps, described elsewhere in the literature [7]. In contrast the upland is much more homogenous and characterized by an intact vegetation layer that consists of a dry moss-, sedge- and dwarf shrub-dominated tundra type and sparsely distributed ponds.

3. Data and Methods

3.1. SAR Data and Processing

The Lena Delta was designated as a long-term monitoring site by the German Aerospace Centre (Deutsches Zentrum für Luft- und Raumfahrt, DLR, Oberpfaffenhofen, Germany) as part of the international collaboration of the PSTG within the World Meteorological Organization. As such, a unique high temporal resolution time series of TSX data acquisitions exists for the central Lena Delta. A time series consisting of 71 SAR images was acquired between 11 March 2013 and 2 October 2015 at approximately 08:34 local time (UTC 22:34), spanning a total observation period of 935 days. The chosen orbit was in descending pass with an 11-day orbital revisit time. The SAR antenna was right-looking and the acquisition incidence angle was 31° in the scene center. This orbit configuration allows for observation of the east facing study site without shadowing effects. The TSX images have an azimuth and range pixel spacing of 0.9 and 2.4 m, respectively, and the scene size is approximately 18 × 56 km. The TSX images were acquired in the Stripmap imaging mode and delivered in processing Level 1B as dual polarized Single-Look Slant Range Complex (SSC) images. Because of a change in the dual-polarization configuration from HH/HV to HH/VV by the DLR in May 2016, we used only the consistent HH (horizontal sent and received) polarization channel for our analyses.

All data were processed using the open source Sentinel Application Platform (SNAP) Version 3.0, provided by the European Space Agency (ESA) [45]. In order to compare the backscatter intensity between individual SAR acquisitions the backscatter reflectivity of each acquisition must be normalized. We calibrated all TSX images and calculated the backscatter intensity (σ°) and converted it to the widely used decibel (dB) scale. We retained the highest detail possible and avoided blurring effects by not performing any kind of image filtering. Consequently, a higher level of noise, which originates from the speckle effect that is present in SAR imagery, remains in the resulting backscatter images. Continuous terrain changes between acquisitions at the cliff top prohibit precise terrain correction using a digital elevation model (DEM) and we therefore used ellipsoid-based terrain correction for all images. The geocoded images remained distorted but were comparable among each other. The expected ground range resolution is approximately 2 m, given the radar geometry of the TSX orbit with an incidence angle of 31° and a slant range resolution of 1.17 in combination with the cliff top topography with slope angles up to 45°. The radar image slant range geometry was transposed to a square image pixel size of 2.4 m and images were geocoded using an ellipsoid correction and nearest neighbor resampling to the Universal Transverse Mercator (UTM) Map projection WGS84 Zone 52 North. Exact co-location of selected images was confirmed using two corner reflectors located on southern Kurungnakh. The TSX images were subset to the extent of the study site (cliff top) to reduce computing time (Figure 2).

Areas of interest (AOI) within the SAR image space (see Figure 1b) were chosen based on field knowledge to analyze backscatter characteristics of cliff and tundra across the study period. AOI ‘tundra’ covered an area of 0.36 km² and was representative of an undisturbed dry moss-, sedge- and dwarf shrub-dominated tundra type. AOI ‘cliff’ consisted of 6 subareas, covering a total of 0.02 km² and were situated on sections of the cliff that showed characteristic high backscatter intensity in the TSX time series. The backscatter intensity value of each pixel within the AOIs was extracted from all 71 TSX images. Because frozen ground conditions and wet snow strongly alter the dielectric properties and resulting backscatter of the cliff surface, we divided the time series into “summer“ (24 acquisitions) and “winter” (47 acquisitions) scenes based on the meteorological dataset. Summer was defined as images acquired when mean daily air temperatures were above 0 °C. Summer images that were likely subject to snow melt or diurnal freezing were not considered for cliff top extraction.

3.2. Automated Cliff-Top Line Extraction from SAR Data

In order to measure the amount of erosion between TSX acquisitions, the cliff tops were automatically generated for all summer TSX acquisitions (Figure 3). We implemented the workflow for cliff top extraction using the ModelBuilder toolbox in ArcGIS™ (ESRI© 10.3.1). The approximately perpendicular and tilted orientation of the active cliff towards the SAR sensor introduces a foreshortening effect to the SAR image that we use for image segmentation. Foreshortening further increases the higher backscatter intensity of the cliff and allows for better differentiation between cliff and tundra. We evaluated thresholds from statistic backscatter values of the AOIs against visually derived thresholds. The threshold of −10.5 dB that resulted from the difference between the means of summer backscatter intensity over the cliff and tundra AOIs did not lead to satisfactory discrimination of the cliff. Ultimately, a range of intensity thresholds was visually evaluated for the segmentation of the upper cliff boundary and tundra and a global threshold of −7.8 dB was chosen as a parameter in the automated cliff top extraction method. Using this threshold, a binary classification was performed on each backscatter image classifying all values above the global threshold as 1 (cliff) and all values below as 0 (not cliff). Due to the remaining speckle effect in the SAR imagery, some pixels within the cliff segment were classified as background and some pixels on the tundra were classified as cliff. To address this problem we applied a two-step neighborhood filtering processing on each binary image based on closing and opening morphological filters. Both were applied using the smallest possible step size of one pixel. The cliff segments of the binary raster files were then converted to polygon vector files. In order to avoid a step-like structure in the cliff-top lines that would result from a strict conversion of the raster cell edges to a polygon, the polygons were simplified based on the algorithm by Douglas and Peucker [46]. We extracted a total of 24 cliff-top lines between 2013 and 2015 from the TSX time series.

3.3. Quantification of Cliff-Top Erosion with the Digital Shoreline Analysis System (DSAS)

In order to quantify cliff-top erosion from the extracted cliff-top lines, we used the Digital Shoreline Analysis System (DSAS, version 4.1) in ArcGIS™ (ESRI©, 10.3.1) provided by the United States Geological Survey [47]. Using DSAS we calculated the distance in m between two cliff-top lines at the intersections with transects that are orientated perpendicular to the cliff section. We used 387 transects with 5 m spacing and a length of 250 m. We extracted cliff-top lines for every 11-day orbital pass but calculated the distance between cliff-top lines at 22-day intervals, rather than the possible 11 days to match the expected erosion at the study site to the TSX ground resolution.

The distance between two time steps is from here on referred to as cliff-top erosion (CTE). CTE was calculated for (1) all 22-day intervals of the TSX cliff-top lines; and (2) the first and last TSX cliff-top line of each year. Negative and positive CTE values reflect retrogressive and progressive cliff-top movement, respectively. Environmental factors have temporal and local differing impacts on each TSX backscatter acquisition and consequently on the derived cliff-top lines. We used the interquartile range of all CTE measurements on one transect within a year to identify outliers. We assumed that the outliers do pick up erosion in most cases but tend to over- or underestimate CTE. If a CTE value was above or below the interquartile range of all CTE measurements along the same transect within a year it was changed and set to the 1st or 3rd quartile, respectively. We assumed only retrogressive movement of the cliff top to be true and set all positive values to zero. Additionally, we set all values that are below 1.2 m, half a pixel size of TSX, to 0, since erosion at this scale is not resolved by TSX (see Section 3.1).

Massive block failures with decameter scale commonly occur at permafrost cliffs with thermo-abrasion through waves [9]. In this case, our filtering could potentially falsely identify cliff top accretion instead of actual erosion if the method detects the cliff-top line of a tundra block that is about to collapse and moves in offshore direction. However, at the study site presented in this paper, eroded tundra blocks are commonly smaller than 1 m and well below the detection limit of TSX, and this effect is therefore not present in the data.

3.4. Cliff Top Mapping from Optical Satellite Data

During the time of the TSX time series we obtained one Worldview-1 scene acquired on 4 August 2014 in panchromatic imaging mode and 0.5 m spatial resolution and one Worldview-2 scene acquired on 26 August 2015 in panchromatic and multispectral imaging mode with a spatial resolution of 0.5 and 1.8 m, respectively (Figure 4). The Worldview-2 multispectral data were pan-sharpened to obtain the highest spatial and spectral resolution using Geomatica software (PCI Geomatic’s© Geomatica 2015 OrthoEngine Service Pack 1). Both images were orthorectified using rational polynomial coefficients (RPC) in an image block adjustment, which uses the geolocation information from both images in the same mathematical model. Using RPCs allows for subpixel geo-positioning accuracies with very high resolution imagery when using only a single ground control point [48]. We used a network of high accuracy ground control points from a field survey in 2013 for precise geolocation with subpixel accuracy in the UTM Map projection WGS84 Zone 52 North. The study site presents a challenging environment for delineation of high-precision ground control points since the tundra environment is very dynamic and does not feature static features such as road intersections or buildings. We used characteristic features such as small ponds and intersections of polygonal ice wedge networks which appeared stable in comparison to older high spatial resolution imagery as ground control points. As topographic input we used a digital elevation model derived from a panchromatic stereo image triplet of the ALOS PRISM Sensor with acquisition dates of 21 September 2006 with an absolute height accuracy of 4.41 m [49]. The cliff top position was manually mapped from both images at a mapping scale of 1:500, which proved to be appropriate for capturing fine details of the cliff-top line [14]. Because the cliff top could not be clearly identified in the northern part of the cliff, the mapping of the optical data was only performed in the area of DSAS Transect IDs 1 to 296 covering a length of 1475 m. DSAS was used to quantify CTE between the two optical cliff-top lines.

3.5. In-Situ Observations of Cliff-Top Erosion

We established a 15 m transect perpendicular to the cliff top to record in situ erosion in 2015 (Figure 5). The transect location was chosen based on site knowledge and optical satellite imagery, which showed sustained rates of erosion and a representative cliff morphology. A total of 29 equally spaced (50 cm) wooden stakes were installed and monitored using a battery powered Brino TLC200 pro time-lapse camera mounted on a 2 m aluminum pole. Images were acquired sequentially every four hours. The camera was set up on 21 June 2015 and operated until 30 August 2015, capturing a total of 70 days. The date of the first visible instability of a stake was recorded as the date of erosion. We then calculated a theoretical daily cliff-top erosion rate.

3.6. Climate Data

Micro-meteorological data records from a long-term measurement station on Samoylov Island (72°22′26′′N and 126°29′45′′E), approximately 8 km northeast of the study area, were used to investigate the influence of climatic variables on erosion (Figure 6). The station belongs to the Fluxnet network with the ID RU-Sam (http://fluxnet.fluxdata.org/). Aggregated 30 min means of air temperature at 2 m height (°C), precipitation (mm water column), and incoming shortwave radiation (W/m²) were resampled to daily means. When less than 75% of the total measurements of the day were recorded the daily mean value was omitted. Thawing-degree days (TDD) of each year were calculated by summing positive daily mean temperature values. To examine the influence of climatic variables and their interactions on the intra-annual cliff-top erosion, we compared the median erosion rates of each of the 22-day TSX intervals from the DSAS analysis (see Section 3.3) with corresponding 22-day aggregated mean values of wind speed, temperature, TDD, incoming solar radiation as well as 22-day aggregated sums of precipitation.

3.7. Statistical Data Analysis

To determine inter-annual variation of CTE and climate data we used analysis of variance (ANOVA). Mean CTE and net CTE of the entire cliff top, as well as mean temperature, TDD, mean incoming shortwave radiation, mean precipitation, precipitation sum, and wind speed were compared by year. Tukey’s Honest Significant Differences (HST) post-hoc test was used to determine where the significant differences were for the annual data. Annual CTE and climate variables were subsequently qualitatively compared. We compared the recorded in situ erosion data to the TSX CTE measurements on the DSAS transect that corresponded to the in situ transect. A simple t-test was used to determine if the in situ and corresponding TSX cliff-top erosion data were significantly different. A significance level (p) < 0.05 was used for all analyses. All analyses were performed using the stats package in the R Studio software environment [50].

We used linear mixed models implemented in the lme4 package within RStudio to examine the intra-annual variation of CTE within each year [51]. The models were built using the 22-day CTE values with corresponding 22-day aggregated climate variables. Linear mixed models were chosen as they permit the regression of multiple explanatory variables (fixed effects) while also accounting for biases associated with time series data such as repeated measures (random effects). We defined wind speed, temperature, incoming solar radiation and precipitation as fixed effects and transect ID and mid-date as random effects to account for repeated measurements. Models were built for each year with increasing complexity beginning with the null model (no fixed effects), to all fixed effects with no interactions, and finally all fixed effects plus interactions. Models were fitted using restricted maximum likelihood due to differing fixed effects. The Akaike Information Criteria (AIC) is a criterion of the quality of the model, with lower AIC values showing better models. An ANOVA was conducted between the best fixed effect model and the null model to determine if the two were significantly different.

4. Results

4.1. TSX Erosion versus In-Situ and Optical Datasets

Extracted TSX cliff-top lines using the backscatter intensity threshold from summer 2013 to 2015 showed gradual and constant erosion at sites with simple cliff geometries (Figure 7a). At locations with complex cliff geometries (gullies) and incised hinterland surfaces strong deviations up to tens of meters from previous or subsequent cliff-top lines were observed. This cliff geometry-dependent performance was also reflected in the comparison with in situ and optical data (Figure 7b). The accumulated erosion of the in situ and two TSX datasets (TSX1 and TSX2), each representing 22-day intervals without 11 day temporal overlap, at a cliff section with simple geometry, were in the same order of magnitude (

Table 1. Comparison of in situ and TerraSAR-X derived cliff-top erosion (CTE) and theoretical cliff-top erosion rates at the location of the in situ setup. Erosion rates were calculated by dividing the net CTE by the number of days with a mean daily temperature above 0 °C in each observation period. The presented values are single measurements.

Sensor	Observation Start (yyyy-mm-dd)	Observation End (yyyy-mm-dd)	Days above 0 °C	Net CTE (m)	Theoretical Cliff-Top Erosion Rate (m/day)
In situ	2015-06-21	2015-08-30	71	−7.0	0.09
TSX	2015-06-25	2015-08-30	67	−5.19	0.07

) and a t-test of theoretical daily erosion rates confirmed that the datasets were not significantly different (p < 0.05).

When TSX and the optical data were compared across the entire cliff top including areas with complex geometry, there were greater deviations across the complex cliff sections (Figure 8). Overall, the southern half of the cliff (Transect ID 1–149) was characterized by stronger erosion and simple cliff geometry while the middle and northern part (Transect ID: 150–300) was characterized by complex geometry and higher variability in erosion. This spatial geometry and variability is reflected in the comparison of both datasets (Figure 8a). Some of the disagreement between optical and TSX data could also be due to the differing methods of terrain correction causing shifts of ±1 transect, equaling ±5 m horizontal offsets. Despite these deviations, the net CTE and theoretical erosion rates of the optical and TSX data were overall consistent (

Table 2. Comparison of optical and TSX-based cliff-top erosion (CTE) measurements. Erosion rates were calculated by dividing the net CTE by the number of days with a mean daily temperature above 0 °C in each observation period. The presented values are single measurements.

Sensor	Observation Start (yyyy-mm-dd)	Observation End (yyyy-mm-dd)	Days above 0 °C	Mean Net CTE (m)	Theoretical Cliff-Top Erosion Rate (m/day)
Optical	2014-08-19	2015-08-26	126	−5.78 ± 4.45	0.04
TSX	2014-08-21	2015-08-30	128	−6.92 ± 6.31	0.05

).

4.2. Inter- and Intra-Annual Cliff-Top Erosion and Climate Data

Inter-annual variability of CTE and climate data is shown in

Table 3. Annual characteristics of erosion and climate at the study site. The mean and median cliff-top erosion (CTE) per 22-day represents the average of all TSX 22-day steps within one year. In contrast, the mean net CTE is derived from the comparison of the first and last TSX summer observation of each year. The climate data in this table represent a temporal subset from June to September.

	2013	2014	2015	Long Term Mean *
Mean CTE per 22-day (m)	−1.37 ± 1.55 a	−1.96 ± 1.85 b	−1.85 ± 1.97 c
Median CTE per 22-day (m)	−1.24	−1.9	−1.49
Mean net CTE (m)	−4.1 ± 3.7 (88 days) a	−6.9 ± 4.6 (88 days) b	−5.1 ± 5.3 (66 days) c
Mean Temperature (°C) (June to September)	6.1 ± 4.8 a	7.6 ± 4.7 b	6.4 ± 4.4 ab	6.4 ± 3.4 (1998–2014)
Thawing-degree days (June to September)	691.9	824.2	628.9	790 ± 120 (1998–2014)
Mean incoming short wave radiation (W/m2) (June to September)	133 ± 87.2 a	141.4 ± 91 a	130.9 ± 81.7 a	-
Mean Precipitation (mm) (June to September)	0.9 ± 2.6 a	0.8 ± 1.6 a	1.1 ± 2.9 a	1.3 ± 3.1 (1998–2008)
Precipitation sum (mm) (June to September)	109	104	146	124 ± 69 (1998–2008)
No of extreme precipitation events >15/10/5 mm (June to September)	2/3/5	0/0/5	2/5/8	1.45 ± 1.30/2.58 ± 1.83/ 6.5 ± 4.23 (1998–2008)

^a,b,c different letters within a row represent significantly different means with p < 0.05; * Long term mean data (1998 to 2012) was taken at a micro-meteorological weather station [34] with the same instrumentation as the climatic record.

. The mean CTE of all 22-day TSX periods was highest in 2014 with −1.96 m ± 1.8 m and significantly different in all years (p < 0.05). In 2015 and 2013, median CTE were not significantly different. In addition to significantly higher CTE in 2014, mean temperature, TDD, and incoming shortwave radiation were all above the long-term average while recorded precipitation was lowest and no extreme precipitation events were recorded. In 2015, mean temperatures were close to the long-term mean but TDD were considerably lower and incoming shortwave radiation was at the lowest of the three years. The year 2015 also had above-average precipitation and a high number of extreme precipitation events compared to 2014 and 2013. In 2013 and 2015, between five and seven precipitation events (>10 mm) were detected, while 2014 had only light precipitation (<5 mm). Incoming shortwave radiation was highest in 2014, but in general did not show large variations across the years.

In general, intra-annual CTE did not show any strong or consistent seasonal trends (Figure 9). In 2013 CTE was lowest and could only be detected in the early and late season. In 2014 two peaks of CTE were observed, the early season showed the highest CTE between end of June and mid-July, after which July showed low CTE values. Early August showed high CTE values and late August again had low values. 2015 showed high CTE in the early and late season, with a minimum observed in July.

The influence of climatic variables on intra-annual CTE was examined by linear mixed models and results are presented in

Table 4. Results of the linear mixed models fitted by restricted maximum likelihood with middate (MD) and TransectID (TID) as random effects. Fixed effects included: CTE: cliff-top erosion, SWIN: incoming shortwave radiation, P: precipitation, TDD: thawing-degree days. AIC: Akaike Information Criterion, BIC: Bayesian Information Criterion. p-values represent results from analysis of variance between the model and null model (no fixed effects).

Year	Fixed Effects	Random Effects	AIC	BIC	p-Value
2013	CTE~SWIN + P + TDD	MD + TID	7138.6	7138.6	<0.01
2014	CTE~P	MD + TID	10,993.4	1102.9	>0.05
2015	CTE~SWIN + P	MD + TID	10,962.7	10,998.0	<0.05

. The results of the linear mixed models suggest that overall the measured climate variables do not have a consistent influence on variations in 22-day CTE within each year. In 2013 the combination of incoming shortwave radiation, TDD, and precipitation sum had the strongest effect on erosion. In 2014, precipitation sum alone was identified as an influential parameter on erosion, although the result was not significantly different from the null model (no fixed effects). In 2015 both precipitation sum and incoming shortwave radiation influenced CTE. In 2013 and 2015 the models were significantly different than the null model, while the 2014 model was not.

4.3. Backscatter Time Series

To understand how the surface properties driving erosion rates were changing, the backscatter intensity of the observation period was examined. Backscatter intensity of the cliff was higher in the summer months and showed a decreasing trend from early to late summer (Figure 10). In all years the cliff showed a peak in backscatter intensity at the beginning of summer between mid-April and mid-May during the transition to above freezing air temperatures and snowmelt. The early season peak in 2013 was less strongly expressed than in 2014 and 2015. Also, in 2013 the earliest onset (early May) of temperatures above 0 °C was observed. By contrast, 2015 temperatures did not reach above 0 °C until the beginning of June. The summer tundra backscatter intensity was slightly lower in 2013 and 2014 than in 2015. In the late season of 2013 and 2014 when temperatures dropped below 0 °C cliff backscatter showed a drop of approximately 3 dB. An ANOVA revealed that backscatter intensity was significantly different between all years for both cliff and tundra with p < 0.005. A general decrease in backscatter intensity of the cliff and an increase in backscatter intensity of the tundra was observed over the three years (Figure 10).

A general divergence of mean backscatter values of the cliff and tundra in summer and convergence in winter was observed. There was a clear difference in the summer acquisitions in all years as shown in Figure 5. The chosen threshold of −7.8 dB was below the median of cliff AOIs and above the 1st quartile of the tundra AOI in all years. The mean summer tundra backscatter of −15.5 ± 6.5 dB was significantly lower and less variable than the mean summer cliff backscatter of −6.4 ± 6.9 dB (p < 0.05) (

Table 5. Mean backscatter over Areas of interest (AOI) by year.

	AOI	2013	2014	2015	All
Mean summer backscatter	Tundra	−16.29 ± 6.53 a	−15.62 ± 6.47 b	−14.81 ± 6.47 c	−15.50 ± 6.51
	cliff	−5.29 ±6.65 a	−6.88 ± 6.93 b	−6.69 ± 7.05 c	−6.39 ± 6.93

^a,b,c different letters within a row represent significantly different means with p < 0.005

). Backscatter of the tundra surface increased from −16.3 ± 6.5 to −15.6 ± −6.5 to −14.8 ± 6.5 dB in the summers of 2013, 2014 and 2015, respectively. The cliff surface showed a slight decrease of backscatter from −5.3 ± 6.6 to −6.7 ± 7 dB from summer 2013 to 2015, respectively.

5. Discussion

Three years of riverbank cliff-top-erosion data from the Lena Delta, Siberia highlight the potential of spatially and temporally high-resolution SAR backscatter time series to monitor rapid tundra disturbance by the degradation of ice-rich permafrost at an inter- and intra-annual scale. This study presents a novel application of TSX X-Band to monitor cliff-top erosion in Arctic tundra.

5.1. Inter-Annual Dynamics of Cliff-Top Erosion

The mean annual net erosion of the entire cliff top derived from the filtered TSX data of −4.1 to −6.9 m per year falls within the range of reported rates of riverbank cliff-top erosion at other Siberian river delta sites (Yana and Indigirka) with comparable cliff height and thermo-denudation-affected systems between 1950 and 1990 [19]. The mean annual TSX erosion rate of −6.9 m in summer 2015 was further validated by the optical data, which showed mean erosion rates of −5.8 m falling within the variance range of one standard deviation (Table 2).

Our analyses showed that mean annual net erosion was significantly different between all years. To better understand the inter-annual differences, we examined mean annual climatic variables in relation to long-term, 10-year means. We found that in 2014 erosion was highest, concurrent with above-average air temperatures, the highest thawing-degree days and shortwave incoming solar radiation, as well as low precipitation. This finding emphasizes the link between air temperatures and ice-rich permafrost degradation and highlights the potential for greater release of material in the fluvial system [52]. Interestingly, in 2015 the erosional activity was significantly higher than in 2013 even though air temperatures were comparable and thawing-degree days and solar radiation were lower. This suggests that the difference between these two years was the result of above-average precipitation both in frequency and intensity. Precipitation has been directly linked to ice-rich permafrost erosion in previous studies [53]. The difference between years and the corresponding climatic conditions suggests a contingent sensitivity of erosion to temperature and precipitation. In other words, precipitation appears to be an important factor when temperatures are at or below average whereas when temperatures are above average, precipitation has less of an influence.

5.2. Intra-Annual Dynamics of Cliff-Top Erosion

Erosion was detected in all analyzed 22-day time steps but no obvious seasonal patterns were observed. The onset of erosion occurs in spring when a threshold of air temperature or ground warming is reached. At this time, water availability from snow melt can lead to increasing erosion potential [54]. Rapid riverbank erosion with high seasonality in spring has been recorded in the lower Lena River [13]. In our riverbank system we could not observe any statistically significant seasonal peaks of erosion that would indicate times with high erosion. The base of the cliff section in our study is not made of ice-rich and highly erodible sediments but is composed of a sand unit that buffers the thermo-erosional activity of the Lena River also during spring flood. This prevents the development of thermo-erosional niches and consequently the rapid erosion by block collapse that is an important component at other Arctic riverbank and coastal cliff systems [11,14,17]. Given that the cliff top is decoupled from the thermo-erosional processes of the Lena River, its erosion likely results from atmospheric forcing and depends on the local sediment and ground ice properties [55], as well as regional climatic conditions [4]. This link explains the relatively gradual thermo-denudation at the cliff top. Interestingly, intra-annual (22-day time steps) erosion was most consistently influenced by precipitation as indicated by the results of the linear mixed models. This result suggests that precipitation has a greater influence over shorter periods of time than temperature. This phenomenon has been documented by Kokelj, et al. [53] who found greater sediment outflow from eroding thaw slumps three to four days following precipitation events. Though precipitation was identified as the most important variable statistically, it is the interaction between temperature, precipitation, and incoming solar radiation that best explains the variability in intra-annual erosion rates.

5.3. Backscatter Dynamics of Tundra and Cliff Surfaces

The backscatter time series clearly demonstrates the ability of SAR X-band data to differentiate between undisturbed tundra and the cliff, justifying the applied thresholding method. In addition, examining the backscatter time series provides insight into changing tundra and cliff surface properties. The observed lower cliff backscatter in 2014 and 2015, compared to 2013, could be the result of a flattened surface geometry and/or partial re-vegetation of the cliff surface. Both effects take place during cliff stabilization of ice-rich permafrost riverbank sites. Kanevskiy et al. [11] showed that rapidly eroding ice-rich vertical riverbank cliffs can stabilize to a flattened cliff with Baydzarakhs (soil remaining after melting of massive ground ice) and vegetation within as little as four years. The ground surface flattens when the cliff top retreats further inland while the cliff base, in our case the underlying sand unit, remains stable. The lower cliff slope angle consequently decreases the effect of foreshortening in the SAR data. Additionally, vegetation is likely to establish relatively quickly on previously eroded and now stable flattened surfaces. Fresh vegetation canopy can increase the fraction of volume scattering and positively affect backscatter intensity [27]. Existing studies on the effect of tundra vegetation on backscatter [23,32,56] do not address the effect of fresh vegetation on previously degraded tundra surfaces. The mixed signals of vegetated and eroding surfaces within the cliff AOIs could explain the higher standard deviations of backscatter in 2015. The active cliff top eroded inland while our chosen cliff AOIs were static and stayed at the same location. Consequently, the scar zone slowly moved out of the AOI and was less represented in the backscatter signature. This time series suggests that segmentation of tundra and cliff using a simple thresholding on X-band backscatter could become more challenging if summer backscatter signatures of cliff and tundra continue to converge over time.

The backscatter time series additionally showed variability during the winter-summer transition, which is very likely driven by changes in the dielectric properties of the ground. When the ground surface was frozen during winter, backscatter signatures of the cliff were lower than in summer but during the transition, values were highly variable over short periods of time. As a result, SAR acquisitions dated before and during snowmelt and during frozen ground conditions were not considered for cliff top retreat analysis. This exclusion could potentially misrepresent important early season dynamics and annual estimates of erosion, but CTE measurements between the last acquisitions of a season to the first acquisition of the following season showed low erosion values, indicating little erosion in this time. In addition, the observed variability of backscatter could indicate the influence of soil moisture on backscatter from bare surfaces in summer which has been shown to increase backscatter by up to 5 dB in X-Band [57]. The impact of soil moisture on backscatter could have implications for image thresholding and resulting cliff-top extraction.

Complex cliff geometries at water tracks and gullies and less steep stabilized cliff sections could also introduce variability into the TSX backscatter data. The sensitivity of radar to non-line-of-sight features is known, terrain features with steep scarps not aligned to orbit direction have been shown to introduce illumination error [58] and shallower slopes impact backscatter angle, both reducing backscatter intensity and the foreshortening effect [59]. These complex non-linear features can introduce errors to the extracted cliff-top line position and should be taken into consideration when monitoring at intra-annual timescales.

6. Conclusions

This study demonstrates that high temporal TSX X-Band backscatter time series in HH polarization can monitor intra-and inter-annual cliff-top erosion dynamics of a Siberian ice-rich riverbank at 22-day intervals. This resolution represents a temporal resolution that was previously challenging for Arctic remote sensing data. A backscatter threshold of −7.8 dB was used to separate the stable and low backscatter signature of −15.5 dB of undisturbed tundra upland and the high backscatter signature of −6.4 dB, of non-vegetated cliff top of summer time acquisitions. We successfully derived cliff-top erosion data across three years from a total of 24 acquisitions. We evaluated these TSX-derived CTE against in situ and optical satellite observations and found that they show comparable seasonal and annual cliff-top erosion. The method is affected by high terrain complexity, by non-line-of-sight features and also by high surface soil moisture. All three impact the automatic extraction of cliff top and make TSX susceptible to over- and underestimation.

We show that ice-rich Yedoma riverbank cliff tops erode consistently throughout the thawing season with little seasonality, largely through thermo-denudation. Maximum erosion rates appeared to be affected by warm air temperatures but precipitation is also important in years when temperatures are below or at long-term averages. Intra-annual erosion rates seemed to be consistently affected by precipitation. However, unmeasured small-scale microclimatic parameters, micro topography and ground ice content are also likely to impact the magnitude of the erosion rates. Overall, this novel method represents a promising new tool to assess rapid changes linked to tundra disturbance and permafrost degradation.