Observing Muostakh disappear : permafrost thaw subsidence and erosion of a ground-ice-rich island in response to arctic summer warming and sea ice reduction

Observations of coastline retreat using contemporary very high resolution satellite and historical aerial imagery were compared to measurements of open water fraction, summer air temperature, and wind. We analysed seasonal and interannual variations of thawing-induced cliff top retreat (thermo-denudation) and marine abrasion (thermoabrasion) on Muostakh Island in the southern central Laptev Sea. Geomorphometric analysis revealed that total ground ice content on Muostakh is made up of equal amounts of intrasedimentary and macro ground ice and sums up to 87 %, rendering the island particularly susceptible to erosion along the coast, resulting in land loss. Based on topographic reference measurements during field campaigns, we generated digital elevation models using stereophotogrammetry, in order to block-adjust and orthorectify aerial photographs from 1951 and GeoEye, QuickBird, WorldView-1, and WorldView-2 imagery from 2010 to 2013 for change detection. Using sea ice concentration data from the Special Sensor Microwave Imager (SSM/I) and air temperature time series from nearby Tiksi, we calculated the seasonal duration available for thermo-abrasion, expressed as open water days, and for thermo-denudation, based on the number of days with positive mean daily temperatures. Seasonal dynamics of cliff top retreat revealed rapid thermo-denudation rates of −10.2± 4.5 ma in mid-summer and thermo-abrasion rates along the coastline of −3.4 ± 2.7 ma on average during the 2010–2013 observation period, currently almost twice as rapid as the mean rate of −1.8 ± 1.3 ma since 1951. Our results showed a close relationship between mean summer air temperature and coastal thermo-erosion rates, in agreement with observations made for various permafrost coastlines different to the East Siberian Ice Complex coasts elsewhere in the Arctic. Seasonality of coastline retreat and interannual variations of environmental factors suggest that an increasing length of thermo-denudation and thermo-abrasion process simultaneity favours greater coastal erosion. Coastal thermoerosion has reduced the island’s area by 0.9 km (24 %) over the past 62 years but shrank its volume by 28× 10 m (40 %), not least because of permafrost thaw subsidence, with the most pronounced with rates of ≥ −11 cma on yedoma uplands near the island’s rapidly eroding northern cape. Recent acceleration in both will halve Muostakh Island’s lifetime to less than a century.


Introduction
Muostakh Island in the southern Laptev Sea is a prominent example (Are, 1988a, b;Romanovskii et al., 2000;Grigoriev et al., 2009) of thousands of kilometres of unstable unlithified coastline along arctic shelf seas (Lantuit et al., 2011a;Overduin et al., 2014).Along this coast, cliffs border marshy coastal tundra lowlands and islands that are underlain by continuous permafrost and composed of continental late Pleistocene ice-rich permafrost sequences called Ice Complex deposits (Schirrmeister et al., 2013).During sum-F.Günther et al.: The disappearing East Siberian Arctic island Muostakh mer, this coast is no longer protected by sea ice and retreats at erosion rates of −2 to −6 m a −1 (Grigoriev et al., 2006).Large areas of up to −3400 m 2 km −1 of coastline are lost annually and currently this rate is more than twice as rapid as historical erosion (Günther et al., 2013a).The distinguishing feature of polar coasts is the presence of a variety of ice types on and ground ice below the earth surface (Forbes and Hansom, 2011).Ogorodov (2011) emphasises the influence of hydrometeorological conditions on the development of coastal thermo-erosion, in particular of thermal and wave energies, both of which are linked to sea ice extent and duration.For the Laptev Sea, Markus et al. (2009) report that the duration of the sea-ice-free season increased on average by 10 days over the last decade, exceeding the average increase of 2 days around the Arctic Ocean.Especially when considering the warming trend of cold continuous permafrost (Romanovsky et al., 2010) and the vulnerability of deep organic carbon to mobilisation (Grigoriev et al., 2004;Grosse et al., 2011), it is important to assess the impact of current climate warming in the northern high latitudes not only as an external disturbance force on ice-bonded permafrost coasts, but also on permafrost-thaw-related land surface lowering of presumably undisturbed adjacent territories.However, such information is practically non-existent for Siberia.
As a consequence of coastal erosion, clastic material enters the near-shore zone (Are, 1998;Jorgenson and Brown, 2005), where it is deposited, reworked and transported (Overduin et al., 2007;Winterfeld et al., 2011).Because ground ice occupies a large proportion of the land's volume above and below sea level, a much smaller amount of material is removed by wave action after thaw than along ice-free coastlines and high rates of coastline retreat are the result (Zhigarev, 1998).Are et al. (2008) conclude that it is mostly thawed material that is being eroded, rather than permafrost.
Coastal thermo-erosion includes two related processes that work temporally and quantitatively differently together.Thermo-denudation (TD) is comprised of the thawing of exposed permafrost, the upslope or inland propagation of a retreating headwall and the transport of material downward to the bottom, all under the influence of insolation and heat flux on the slope (Mudrov, 2007).Thermo-abrasion (TA), on the other hand, is defined as the combined action of mechanical and thermal energy of sea water at water level (Are, 1988a).Despite temporal variations in their intensity, both processes are interconnected, since TD sooner or later becomes inactive after TA comes to a standstill.
Multitemporal applications of remote sensing data are of particular interest for assessing permafrost-related natural hazards such as erosion of frozen sea coasts and thaw subsidence (Kääb, 2008).Numerous recent change detection studies exist and aim to identify coastline variations in different permafrost settings (Lantuit et al., 2013).In concert with time-lapse photography, Jones et al. (2009a) analyse the coastal erosion development around Cape Halkett using high-resolution remote sensing data of the northern Alaska sea coast.They find that, after increasing slightly over the last 5 decades, annual erosion accelerated abruptly and almost doubled, reaching −13.8 m a −1 from 2007 to 2009.They attribute this increase to more frequent block failure as a consequence of higher sea surface temperatures and longer fetch, which potentially create more erosionally effective storm events (Jones et al., 2009b).Lantuit et al. (2011b) study storm climatology and use a set of aerial photographs and satellite images to investigate erosion rates around the entire Bykovsky Peninsula near Tiksi in the Laptev Sea over six consecutive time periods.They show a clear dependency of coastal erosion on backshore thermokarst geomorphology, but do not find either a pronounced temporal trend in the mean annual coastal retreat rate over 55 years (−0.59 m a −1 ), nor a relation to storm activity.For the western coast of the Yamal Peninsula, where retreat rates range from −0.8 to −2 m a −1 , Vasiliev et al. (2006) rely on long-term observational data of the polar station Marre Sale, where the length of the warm period is 102-137 days long, while the open water season lasts for 70 days, on average, generating different preconditions for TD and TA.Although the Kara Sea region experiences frequent storms of long duration (Atkinson, 2005), Vasiliev (2003) finds that only in occasional cases up to 20 % of coastal retreat can be attributed to storms.Arp et al. (2010) report on recent erosion for the Alaskan Beaufort Sea coast, where they observe even more rapid rates of up to −17.1 m a −1 , but find little correlation to sea surface and soil temperatures and, in particular, no consistency with storm events.Although potential local controls on erosion such as ground ice content have been identified (e.g.Dallimore et al., 1996;Vasiliev, 2003), it is difficult to establish a relationship of erosion of permafrost coasts to one or another external factor.Moreover, since current environmental changes are expected to intensify coastal erosion, there is a sustained need for information on coastline recession rates in conjunction with seasonal observations in order to better understand the mechanisms driving thermoerosion and subsequent land loss along permafrost-affected coasts.In addition, however, the influence of coastal thermoerosion on permafrost degradation processes in backshore areas has received scant investigation.
The main objective of this paper is to systematically analyse seasonal thermo-erosion dynamics and backshore degradation for a ground-ice-rich permafrost coast in the Laptev Sea.We use a set of contemporary very high resolution satellite imagery, repeated geodetic surveys in the field and historical aerial photographs to provide current (2010)(2011)(2012)(2013) and historical (since 1951) quantification of planimetric land loss, volumetric coastal erosion, and land surface lowering due to thaw subsidence in backshore areas.In conjunction with digital elevation models (DEMs), we use a geomorphometric method for assessing macro ground ice content of Ice Complex deposits, in order to consider this factor for the estimation of the mass of material that must be reworked by coastal thermo-erosion following thaw and the resulting sed- iment supply to the nearshore zone.Using time series of local sea ice concentration and air temperatures, we apply normalisation to coastal retreat observations over seasonal and interannual periods to identify their seasonal intensity and to discuss environmental controls on processes involved in coastal thermo-erosion development.

Study site
Muostakh is a small island (70 • 35 N, 130 • 0 E), in the Buor Khaya Gulf of the southern central Laptev Sea (Fig. 1), located 40 km east of the harbour town Tiksi in northern Yakutia (Russian Federation).Though situated on the ocean, the severe subpolar climate with mean annual air temperatures in Tiksi of −12.9 • C (1933-2013), where the warmest month does not exceed 10 • C, is continental due to prolonged sea ice cover.Muostakh lies within the northern tundra zone.The vegetation cover is characterised by moss-grass, lichens and dwarf shrub tundra.Cryogenic micro relief features are widespread and include mud boils, frost cracks, peat mounds, thermo-erosional gullies, high-centred polygons on inclined surfaces and thermokarst mounds (baydzharakhs) on coastal bluffs.The island has an elongated narrow form oriented SSE-NNW and is approximately 7.5 km in length with a maximum width of ≤ 500 m at sea level.At the southern margin, next to the former polar station Muostakh, Ostrov, a lighthouse marks the navigable channel into the sheltered Tiksi Bay.As a continuation of the island, an interrupted sand spit chain extends another 5.2 km southwards.Grigoriev (1993) supposed that Muostakh Island was formerly connected with the Bykovsky Peninsula further in the north (Grosse et al., 2007), but nowadays they are separated by a distance of 15.8 km.Both Bykovsky and Muostakh consist of Ice Complex deposits and their sedimentological and cryolithological structures suggest simultaneous formation (Slagoda, 2004).According to the Mamontovy Khayata section on Bykovsky, the Ice Complex in this area formed from 58.4 to about 12.2 ka BP (before present) (Schirrmeister et al., 2002), the clastic material is of local origin (Siegert et al., 2000), and accumulated during the subaerial exposure of the East Siberian shelf.Subsequent Ice Complex degradation through thermokarst resulted in alternating relief of depressions (alas) and uplands (yedoma) (Morgenstern et al., 2011, and references therein).Peat and wood on the base of the Holocene cover on Muostakh showed ages in the range of 2-7 ka BP.Muostakh represents a remnant of the late Pleistocene accumulation plain that remained after the sea level drew to near the current level 8 ka ago (Gavrilov et al., 2006), and the highstand of the Holocene transgression was reached 5 cal.ka BP (Bauch et al., 2001).It serves as a witness for the widespread occurrence of Ice Complex islands on the shelf that have been completely destroyed by coastal thermo-erosion (Gavrilov et al., 2003).Ice-poor sands of Pliocene-early Pleistocene age underlay Ice Complex deposits (Slagoda, 2004).Ice complex thickness on Muostakh is 31 m, 10 m of which extend below sea level (Kunitsky, 1989), providing very favourable conditions for TA.
According to Kunitsky (1989) the permafrost temperature of the non-degraded yedoma on Muostakh at the depth of zero amplitude is −10.4 • C, which is cold permafrost and a typical value for yedoma uplands at this latitude in northern Yakutia (Romanovsky et al., 2010).Since Ice Complex deposition took place under permafrost temperatures of −25 to −28 • C (Konishchev, 2002) (2010), alas temperatures at this latitude are around −9 • C, suggesting that permafrost in the coastal zone has undergone additional thermal degradation.Also in the south of the island, around the former polar station, Slagoda (2004) identified fragments of an alas.Regular stationary monitoring of coastal erosion in the Laptev Sea is conducted only in two places: Mamontovy Khayata on the Bykovsky Peninsula and on Muostakh Island (Grigoriev, 2008).Based on these time series, Muostakh is famous for very high erosion rates, where the northern end of the island for example retreated by about 25 m in 2005 and the nearby east-facing coast by 11 m (Fig. 2).Along with rapid erosion rates, the morphology of the coastal cliff may substantially change its appearance (Fig. 3).

Field work
During a joint Russian-German expedition to Muostakh in August 2011, a network of well-distributed geodetic anchor points was established as a precondition for consistent repeat topographic surveys and their transformation to an absolute coordinate system (Günther et al., 2013b).During a subsequent expedition in August 2012 (Opel, 2015), a repeat survey was conducted.We used a ZEISS ELTA C30 tacheometer for distance and height measurements.Concentrated mainly along the coast, 2392 points were measured that cover about two-thirds of the island's 15 km coastline perimeter.In local project coordinates the point cloud was highly self-consistent, while the absolute geocoding accuracy had a root mean squared error (RMSE) of 1.36 m.Dur-  ing a 1-day visit in 2013, selected points were resurveyed within the long-term monitoring.
Measurements of active-layer thickness (ALT) were conducted during the 2012 expedition for the period from 15 to 23 August.In order to capture spatial variability of ALT, altogether 323 ALT measurements were made by mechanical probing.Mapping of ALT across the island was done along transects of 9 km length with an equidistancy of ≤30 m and comprised all soil and vegetation associations (Fig. 4).

Data fusion and change detection
Remote sensing data were acquired on different dates in order to create a time series of images that was integrated into a geographic information system (GIS), to detect and measure land loss resulting from coastline position changes.In this study, GIS serves as a basis for combining field survey data, historical aerial imagery, contemporary satellite images, and products generated from these data such as DEMs, orthoimages, and digitisation records.
The use of satellite images from different sensors with varying spatial, spectral, and radiometric properties repre-sents challenges for change detection.This is especially true when using very high resolution image data with a ground resolution of < 1 m (Dowman et al., 2012), not only because pixel-based approaches have been designed for lowto medium-resolution imagery (Hussain et al., 2013) but, specifically, because of the different acquisition geometries that must be considered, requiring careful geometric rectification and topographic correction.Although our study area is flat tundra lowland, our main object of interest, the upper coastline, is always located on the sharp edge of coastal cliffs where abrupt elevation changes occur.In addition to temporally very dynamic elevation changes, steep thermoabrasional cliffs and baydzharakhs on thermo-denudational coastal bluffs cause large dynamic shadow effects due to different illumination angles of the low solar elevation in high latitudes.Varying conditions of sea-ice-covered and sea-icefree coastal waters, as well as the presence or absence of banks of snow at the cliff bottom, lead to large reflectance variabilities between acquisitions, making radiometric calibration almost impossible.These conditions lead to problems with automated change detection techniques, examples of which are given in Kääb et al. (2005).
Satellite images must be georeferenced for spatial calibration of multitemporal and multisensor data for change detection.While georeferencing corrects for most distortions connected with the acquisition system, orthorectification corrects for relief-induced displacement effects and creates calibrated satellite image products with the geometry of a map, allowing for distance and area measurements.
GIS-related work was done using ESRI ArcGIS 10.1.Cliff top line positions were manually mapped in orthorectified imagery at different points in time.Unlike other studies on coastal thermo-erosion using the transect method (Günther et al., 2013a), we derived vector data of areal land loss and subsequently calculated seasonal variations of TD.This approach is also used by Aguirre et al. (2008) and Tweedie et al. (2012) for monitoring changes along an arctic coastline but on the basis of DGPS (differential GPS) measurements.Since we also aim to analyse very short time periods of a few days, for which erosion is expected to be rather chaotic, this approach ensures that every event is captured.We divided the studied coastline into 118 segments of 50 m width, with 29 located along the western coast and 89 on the eastern coast.Normalisation of eroded area by baseline length of each segment provided absolute linearised coastal retreat in metres and normalisation over time rates in metres per year.Cliff bottom line position changes are regarded as baselines for historical and subdecadal TA dynamics.We use TD and TA to refer to the rates of coastline position change per year.

Aerotriangulation of historical air photo strip
Aerotriangulation, or block adjustment of a bundle of rays from object to image coordinates, is a standard method in photogrammetry (Konecny and Lehmann, 1984).Aerial pho-

F. Günther et al.: The disappearing East Siberian Arctic island Muostakh
togrammetry is well-suited to quantify historical decadalscale temporal change (Kääb, 2008).Although coastal erosion studies in the East Siberian Arctic using early aerial photography are common (Grigoriev, 1993;Are, 1999;Lantuit et al., 2011b;Pizhankova, 2011), few use stereophotogrammetry, generally because image parameters for old aerial photos are unknown.Therefore, valuable elevation information available from these data sets still remains untapped.
Five airborne images covering Muostakh Island taken on 9 September 1951 along one flight strip were utilised in this study.Hard copies of 180 × 180 mm edge length were scanned using a photogrammetric scanner at 14 µm scan resolution, corresponding to ≈ 0.4 m on the ground.No information on focal length, principal point offset and radial lens distortion was available.However, the latter two can be compensated to some extent with exterior orientation (spatial location of the projection centre and camera's view direction), for which camera focal length is required (Jacobsen, 2001).Whether or not the correct focal length of the air survey camera is used, calculation of the flight altitude is necessary (Knizhnikov et al., 2004), which can be done by determining the scale of the frame photography.The scale number S a of the frame photography was roughly estimated following by measuring the same distance between two objects within the original photograph (d a , cm) and on-screen (d m , cm) within a contemporary orthoimage with a map scale (S m ) set to 1 : 10 000.This approximation resulted in a large scale of 1 : 28 000.Using the Aerial Photography model in PCI Geomatica's 2013 module OrthoEngine, we collected fiducial marks in each image to visually define the principal point, collected a set of 25 stereo GCPs (ground control points) and automatically computed over 850 tie points (TPs) using cross-correlation for strip stabilisation.Bundle block adjustment was performed iteratively with the focal lengths of air survey cameras that existed at that time according to Shcherbakov (1979).The best overall solution was achieved with a focal length of 100 mm (likely Liar-6, 104 • wide angle lens, used for topographic medium-scale mapping), yielding a RMSE for GCP locations of 2.4 m and TPs of 0.3 m.For the exterior orientation parameters, this corresponds to a flight altitude of around 2600 m, which is consistent with the theoretical flight height of 2800 m calculated from the approximate photo scale following where c k is the focal length and h flight height, both in metres.The air survey strip constellation provides alongtrack stereo and triple overlap situations, with base-to-heightratios of 0.7 and 1.4, respectively, allowing for height parallax measurements and DEM extraction of steeper slopes as well as over flat terrain.Vertical accuracy of the DEM was estimated to be 2.3 m, based on the difference of input and calculated elevations at stereo GCP locations.Rectangular corners of footprints verify nadir viewing geometry and a robust model, which did not result in overfitting outside the GCP cloud, an important precondition for reconstructing the former shape of the island and subsequent change detection.
Based on the DEM of 1951, the aerial photographs were orthorectified and stitched together in a seamless orthomosaic.

Multisensor block adjustment
Multisensor data fusion, in this study of GeoEye (GE), QuickBird (QB), WorldView-1 (WV-1) and WorldView-2 (WV-2), offers the opportunity to merge images collected from different satellites and different orbits in one triangulation process.According to Toutin (2004), the simultaneous solution of an entire image block offers several advantages, for example the number of GCPs can be reduced, better relative accuracy between images can be obtained and finally more homogeneous orthoimages over large areas can be produced.Satellite sensor models described by rational polynomial coefficients (RPC) provide a high potential of simple and accurate geopositioning (Fraser et al., 2006), are ideally suited for block adjustment of narrow field of view sensors (Grodecki and Dial, 2003), but require some bias correction (Fraser and Ravanbakhsh, 2009), and generally serve only as an approximation of physical sensor models when orbital information is not provided in the metadata (Poli and Toutin, 2012).We performed block adjustment and subsequent orthorectification using our own ground control within the Rational Functions model (RPC-based) in OrthoEngine.
Prior to further geometric processing, we applied singlesensor and single-date image fusion to QB, GE, and WV-2 imagery using the enhanced pan-sharpening method of Zhang (2004).Nine very high resolution images were acquired as standard/orthoready products (Table 1), with panchromatic (PAN) imagery resampled with sinusoidal kernels, for better representation of sharp features (Toutin, 2011).Due to varying moisture and illumination conditions between acquisitions, we found that not all GCPs collected in the field could be identified unambiguously in each image, resulting on average in about four GCPs per image.The RPC model is a viable alternative for rigorous sensor models (Cheng et al., 2003), and several studies show that the effect of the number of GCPs on 3-D RPC block adjustment is limited, yielding almost no further improvement, if configurations of more than four GCPs are used (Fraser and Ravanbakhsh, 2009;Aguilar et al., 2012).Based on our topographic reference measurements, 37 elevation TPs were additionally incorporated into the block, to achieve higher redundancy in RPC bias correction and mainly to better align images to each other.The zero-order polynomial turned out to be the most stable and best possible solution of the block, yielding a submetre accuracy within the entire block of 0.81 m RMSE (Table 1).According to Günther et al. (2013a), the mutual RMSE of each data set (Table 1) is then considered to be the relative georeferencing uncertainty in the determination of the cumulative uncertainty in coastline position, which results out of the combination of this error, the ground resolution, and the additional 2-D positional error introduced by the DEM used for orthorectification, depending on the incidence angle.Uncertainties in change rate calculation over six periods between 15 June 2010 and 17 July 2013 are also applied following Günther et al. (2013a) and were in the range of 0.47-0.65 m.

True orthorectification
The 2013 GE images were acquired as a stereo pair and used for DEM extraction using the stereo model of the entire image block.We applied a DEM editing procedure in order to remove noise, interpolate areas of unsuccessful matching, and for low-pass filtering.The final DEM features high detail and its spatial resolution is 1 m.Using 1158 survey points from 2011 that were located within the 2013 cliff top line extent, the mean elevation difference was 0.06 m, which reflects the good match of both data sets in absolute reference height.Accordingly, the vertical accuracy of the 2013 DEM was evaluated as σ = 0.64 m.However, the 2013 DEM could not be used for orthorectification of the 2010-2012 satellite images, because of its mismatch in cliff top line position compared to earlier dates.
Due to this, radiometric similarities for multidate singlesensor (e.g.QB-QB) and multidate multisensor (e.g.WV-GE) constellations were evaluated with regard to stereoscopic interpretation at an earlier date.Through pansharpening of GE's near-infrared band (0.45-0.8 µm) with the PAN band (0.78-0.92 µm wavelength), the spectral range was adjusted to WV-1 PAN imagery (0.38-0.88 µm) to achieve enhanced image matching.Epipolar image matching of GE acquired on 13 July 2010 and WV-1 acquired 26 days later resulted in visually good topography and showed the least standard deviation between input and calculated elevations for 3-D TPs.
Final DEM generation was performed on a hybrid vectorraster data basis.Contour lines from the stereoscopic DEM, cliff bottom and top lines (digitised in ellipsoid-based orthoimages considering geoid height offset), and all point data of the 2011 and 2012 topographic field surveys (that are necessarily within the 2010 cliff top line extent), were incorporated into a terrain interpolation procedure according to Hutchinson and Gallant (2000).The final DEM represents the island's state in the early summer of 2010.The vertical accuracy of the DEM at survey point locations was σ = 0.3 m, therefore introducing ≤0.1 m 2-D positional uncertainty along the cliff top line through subsequent orthorectification, instead of ≤7.2 m random terrain-induced displacement for each of the very high resolution images when using the initial RPC reference height of the image products.

Elevation difference uncertainty assessment
Differencing of multitemporal DEMs was done using the 2013 and 1951 data sets.The relative vertical uncertainty between both DEMs in the island's interior was −0.4 ± 2.2 m.However, according to the strategies of Nuth and Kääb (2011) and based on our reference data from topographical surveys, a height-dependent bias in DEM difference of 0.47 m m −1 could be identified.Günther et al. (2012) also report systematic underestimations of height measurements in historical aerial photography stereo pairs.DEM errors may result from inexact matching and a lack of contrast within and similarity between stereo images (Nuth and Kääb, 2011).The systematic bias was corrected using an empirical equation of second polynomial order (Fig. 5).
The vertical accuracy of the 2013 and 1951 DEMs was determined to be 0.64 and 2.3 m, respectively.Regarding the 1951 DEM, this is not particularly meaningful, because reference heights were derived from the survey data, collected 60 years later.Restricting the search matrix size for crosscorrelation of the 1951 aerial photographs to a modern elwww.the-cryosphere.net/9/151/2015/The Cryosphere, 9, 151-178, 2015 evation maximum of 21 m a.s.l.resulted in an almost complete failure of height parallax measurements north of the alas.For this reason, we did not restrict the search matrix size and obtained a consistent 1951 DEM with maximum elevations of 24.9 m a.s.l., which is in accordance with a topographic map from 1953.The observed maximum elevation change is −25.5 m, which illustrates that bias correction is robust also in former island areas not covered by contemporary reference data.Since almost the whole terrain of Muostakh may potentially be susceptible to permafrost thaw subsidence or other causes of elevation change, no long-term reference height points were available, except for two elevation indications of 20.9 and 25 m a.s.l. in the 1953 map.The positional inaccuracy of 75 m related to the topographic map was evaluated using the centroids of four lakes that were visible in the 1951 imagery and marked in the topographic map, but gradually drained during the 1960s and 1970s.With respect to the flat relief on yedoma uplands, elevation mismatch was evaluated within buffers of 150 m diameter around both points and revealed an absolute DEM difference uncertainty of −1.56 ± 0.78 m.Within the already degraded alas depression at 3-4 m a.s.l., where no substantial changes are expected, elevation differences were −1.01 ± 0.55 m on average.Together, both indicators suggest a small residual height dependent bias of −0.028 m m −1 and revealed a mean error of 1.31 m that possibly corresponds to a slight overestimation of negative terrain height changes.

Environmental parameters
We use environmental observations to relate coastal dynamics to its potential drivers of atmospheric warming and sea ice reduction.As a proxy for marine abrasion along the cliff bottom line, we use the inverse of daily sea ice concentration that is open water extent per day.To relate the rate of thaw along the cliff top line, we use air temperature (T air ) data as positive mean daily temperatures.These data are then used for correction of coastal erosion rates from different image acquisition periods.To identify single wind events and general wind patterns during the period of open water that potentially influenced coastal retreat we use wind data.

Sea ice concentration data
Daily percent sea ice concentrations from 1992 to 2013, based on Special Sensor Microwave Imager (SSM/I), were used to calculate the open water fraction in percentage per day.Derivation of total sea ice concentration from SSM/I data uses dual polarisation measurements with the 19 and 37 GHz channels of 25 km spatial resolution since 1978.Higher resolution of 12.5 km is available using the 85.5 GHz high frequency channel that did not work before 1992.Although the 85.5 GHz product covers a shorter period of time and might be affected by larger uncertainties over lower sea ice concentrations and open water (Lomax et al., 1995), 25 km resolution is too coarse to study the coastal zone around Muostakh.We worked with the 12.5 km data product of the ARTIST (Arctic Radiation and Turbulence Interaction STudy) sea ice algorithm (Kaleschke et al., 2001) distributed by Ifremer/CERSAT (2000), which is based on a hybrid model and provides reliable results (Ezraty et al., 2007).
In very few cases, the data can also accept negative and positive values outside the 0-100 % range (Andersen et al., 2007), requiring correction.Data was masked by land mass and limited to a 100 km radius around Muostakh Island, spatially corresponding to the Buor Khaya Gulf (Fig. 6).Daily sea ice coverage was smoothed using a 7-day running mean and converted to data of open water fraction in order to determine and count open water days (OWDs).

Air temperature data
The spring to fall seasonal cycle in the Lena Delta region features rising T air during spring until the end of snowmelt, T air well above the freezing point during summer, and fall is characterised by the beginning of refreeze (Langer et al., 2011).In order to evaluate the response of TD to T air over time, we use positive mean daily T air and positive degree-day (PDD) sums (Braithwaite, 1995), obtained from the temperature curve integral above 0 • C in Kelvin days (Kd) (Jonsell et al., 2013).The hydrometeorological observatory Polyarka near Tiksi (WMO # 21824) measured T air 3 times a day from 1932 until 1936, every 6 h until 1970, and has measured it since then every 3 h.Data were downloaded from the electronic archive of observations at Tiksi (Ivanov et al., 2009a, b).For different erosion observation time ranges within the period from 1951 to 2013, PDD were calculated using mean daily T air data as annual sums and as PDD sums over certain seasonal and interannual periods.

Wind data
Wind blowing over open water generates waves breaking on the shore face.Wave height, energy and erosion potential of TA is proportional to wind speed.Wind data of Tiksi from 1992 until 2012 were examined for direction, speed, and single strong wind events.We focus on this period to relate wind data to daily sea ice conditions.Based on measurements every 6 h, mean wind speeds were classified into eight wind directions and analysed for observations during sea ice breakup and the sea-ice-free period of a particular year.

Local parameters
Permafrost is prone to thawing because its core element is the occurrence of ground ice.Ice wedges constitute a large fraction of the subsurface volume.They extend in different generations and stratigraphic units from the top of the permafrost to below sea level, and their size directly determines how much clastic material must be thawed and subsequently removed by coastal thermo-erosion.On Muostakh, they are syngenetic ice wedges.Visually estimating macro ground ice content based on the fractional ice-wedge volume (for example, from photographs of the coastal cliff) is often complicated by slope and perspective, by debris material that obscures undisturbed in situ material, as well as by the fact that the collapse of thermo-abrasional cliffs occurs along ice- wedge axes (Are, 1988b), implying that an exposed wall of ice may not serve as a representative random test for the geological subsurface.

Macro ground ice
Baydzharakhs are a characteristic ephemeral cone-shaped thermokarst landform and represent the remnant frozen sediment core and geometric centres of thawed ice-wedge polygons (Mudrov, 2007).During field work in 2011, twenty coastal slope profiles were surveyed between and across baydzharakhs at different locations, where we observed that baydzharakh spacing varied and might do so depending on their fractional volume of the subsurface (Fig. 7).
Using terrain-corretced satellite imagery, we extended baydzharakh mapping to erosional coast segments (Fig. 8).Taking baydzharakh centres as seeds, we subdivided the surface into cells of a Voronoi diagram (Reem, 2010), which we use as an estimate of polygon morphology prior to thaw.The largest possible circle within a cell was calculated using the maximum Euclidean distance of each cell as radius.Based on the assumption that the sediment centre of each polygon has a cylindric form, macro ground ice content as fractional volume (V wm ) is calculated from where A circle is the surface of the sediment centre and A polygon the total polygon area in square metres.

Subsidence potential
According to Katasonov (2009), the porosity of Ice Complex is very large, due to excess ice and fine particle size.However, natural sediment deposition forms cavities, and this fraction of the total porosity must be disregarded in terms of subsidence.We assume a porosity of 0.4, following where ρ is mean bulk density of the Ice Complex on Muostakh of 1.6 ± 0.25 × 10 3 kg m −3 , according to Solomatin (1965), and ρ O the particle density of non-porous clastic material (2.65 × 10 3 kg m −3 , Strauss, personal communication, 2012).Assuming that all pores are filled with ice (Strauss et al., 2012), the pore ice fraction of intrasedimentary ice does not contribute to subsidence.Therefore, the relative subsidence potential of thawing Ice Complex deposits was calculated following modified after Mackay (1966) and Are (2012), where V wm is volumetric macro ground ice, V s volume of the sediment part, W is intrasedimentary ground ice and φ porosity.

Historical erosion development
The development of thermo-erosion on Muostakh and its shaping of the island was analysed over a period of time of more than half a century.Starting with the first aerial photographs from 1951 and ending with the most recent GeoEye image of 2013, erosion was quantified.We analysed areal land loss and the associated volumetric land losses over 62 years for the entire island.Based on this, we concentrated on the eroding portion of the coastline, where 118 coastline segments were studied in more detail.Each segment corresponds to a 50 m coastline length at beach level.Squares used for symbolising erosion in Fig. 9 (middle) show the spatial distribution of all studied coastline segments.

Mass movements
The volume of Muostakh Island has decreased by 40 % between 1951 and 2013, based on multitemporal DEMs calculated for those years.According to the overall volume change (Table 2), the calculated rate of mean annually eroded volume on Muostakh is 0.45 × 10 6 m 3 a −1 , corresponding to 0.36 × 10 6 t a −1 of annual ground ice thaw and 0.16 × 10 6 t a −1 of sediment displacement.About two-thirds of this results from erosion on all sides of the island and about one third from an overall degradation through surface lowering.
The cumulative eroded volume within the coastal segments studied in detail was 16.3 × 10 6 m 3 , which means that 90 % of coastal erosion was recorded by our coastline subsample (Fig. 9, middle).Within each of the the 50 m coastal segments, annually eroded volumes per segment varied broadly and were on average −2240 m 3 a −1 .Very high values of up to −28 300 m 3 a −1 were detected at the north cape.Generally, eroded volumes were close to the median of −1280 m 3 a −1 , because most of the coastline is eroding more slowly and has a lower backshore height than the northern cape.

Land subsidence and active layer thickness
In addition to large elevation decreases along the eroded coastline, we observed land subsidence across the entire island (Fig. 9, right).The mean elevation of Muostakh in 2013 was 14.4 m a.s.l.The 1951-2013 DEM difference raster was clipped to the interior of the cliff top area of 2013, in order to exclude the influence of coastal cliffs for further analyses of this phenomenon.The mean elevation change was −3.56 ± 1.8 m.Based on these data, except for a very limited area around the former polar station, the island experienced land subsidence at a mean rate of −5.8 ±2.9 cm a −1 .In particular, in the northern part elevation decreases were large  and land subsided by around −10.9 ±0.6 cm a −1 over the last 62 years.Quite unexpectedly, the spatial pattern shows land subsidence is more active close to erosive parts of the coastline and intensifies the more rapidly coastal thermo-erosion is proceeding (Fig. 9, right).This becomes particularly evident not only in the north where rapid rates of coastal erosion coincide with strong subsidence, but also in the middle part of the island, where erosion from both sides leads to stronger subsidence compared to neighbouring non-erosive coastal segments and the adjacent hinterland.
Based on 326 equally spaced transect measurements, ALT on Muostakh in 2012 was on average 47 ± 19 cm.Generally, large ALTs were clustered on well-drained slopes close to the coast, on slopes of the alas, and on the southern tip of the island (Fig. 4), the only place not affected by subsidence.Land subsidence observations over the historical time period were then linked to current measurements of ALT (Fig 10).In contrast to the predominant image of larger ALT causing permafrost thaw and subsequent subsidence, we found, for example, that a shallow ALT of ≤20 cm is associated with intensive subsidence of around −7.5 ± 0.4 cm a −1 , while at locations of deep ALT ≥ 80 cm mean subsidence was only −3.4 ± 1.9 cm a −1 .

Coastline changes
TA was analysed over 61 years, where the start and end points of the observation period are in early September, meaning there is no shift with respect to season.Günther et al. (2013a) found that TD and TA along Ice Complex coasts are interconnected: TA is the limiting component for coastal thermo-erosion intensity on the longterm scale.TA also better reflects the overall land loss of the base area of the island.The base area extent of Muostakh in 1951 was 3.8 km 2 .By 2012 it had shrunk by −23.7 % to 2.9 km 2 , which corresponds to a mean land loss of −14 700 m 2 a −1 .Given Muostakh perimeter of 15.5 km, this corresponds to a mean coastline retreat rate due to TA of −0.95 m a −1 when examined for the entire island including non-erosive coastline sections.
For erosive sections, areal land loss was mapped within 118 segments.Transformation of 2-D areal data to mean distance measurements was done individually via the baseline length of a particular segment.The uncertainty of cliff bottom position change is ±1.7 m, for TA ±0.04 m a −1 .Along the 5.8 km coastline covered by our segmentation (corresponding to the former 6.4 km in 1951), absolute TA had a mean of −109.7 ± 80.6 m, whereas annual rates were in the range from −0.2 to −7.2 m a −1 with a −1.8 ± 1.3 m a −1 mean.According to Fig. 11 (left), 71 % of TA rates over the historical period were clustered towards slower rates of ≥ −2 m a −1 ; for comparison, only 59 % of modern TA rates were slower than the mean value, suggesting coastal erosion over the historical period progressed relatively uniformly, or that temporal averaging of several erosion events occurred.Although the northern cape is eroding at a different angle than the rest of the north-eastern coastline, it has traditionally been of interest as it reflects the increasing distance between Cape Muostakh on the adjacent Bykovsky Peninsula and Muostakh Island.Tracing the position change of the exposed northern cape of the island between 1951 and 2012, maximum absolute cliff bottom line recession was −585 m, which is equivalent to −9.6 m a −1 .

Interannual and seasonal erosion development
Open water and positive T air are unequal in duration and time of year (Fig. 12).In 2010, the first two images were acquired during ice melt and increasing mean daily T air (15, 29 June).Sums of PDD and OWD were correlated (Fig. 13) and are probably generally correlated, since sea ice melt is driven to a great degree by heat exchange with the atmosphere.This means that we expect TD and TA to be correlated insofar as they are driven by PDD and OWD, respectively.Due to the strong seasonal constraints on the development of coastal thermo-erosion, the discrepancy between the start and end points of the observation periods and the duration of the season when TD and TA are able to proceed may result in an over-or underestimation of rates.In cases of mismatch between two acquisition dates, instead of direct change rate calculation only, we corrected calculated rates over time using a season factor.Season factors were derived from the ratio of either the number of days of open water or of positive mean T air during the specific observation period to a perennial reference period.Season factors are used to calculate the actual coastal erosion velocity over a particular period of time and to compare velocities between periods.

Current thermo-abrasion
We examined current dynamics of TA using GE images of 13 July 2010 and 17 July 2013 as the data set spanning the longest period of the recent past, for which the cliff bottom is free of snow.The base area reduction of the entire island was −22 300 m 2 a −1 during the last 3 years (compared to −14 700 m 2 a −1 on average from 1951 to 2012; Table 3).Of this erosion, 89 % occurred within our 118 coastline segments for detailed study.During the last 3 years, mean TA was −3.4 ±2.7 m a −1 and therefore currently 1.9 times faster than over the historical time period (Fig. 11, left).Of all segments, 19 % experienced slight deceleration, while only at a few segments TA rates remained almost unchanged.However, of note is the fact that this almost doubling of coastal erosion is not due to outliers, but derives from a broad acceleration at segments previously eroding in a narrow range from −0.5 to −2.5 m a −1 to currently −1.5 to −8 m a −1 .Between July 2010 and 2013 the northern cape retreated at −51 m, which corresponds to a TA rate of −17 m a −1 (compared to −9.6 m a −1 over the historical period).Erosion at the cape determines the dynamics of the sand spit forma-tion next to it.At this location, additional topographic survey data of the expeditions in 2011 and 2012 were available and covered down to the cliff bottom (Fig. 14).Interannual variations in coastline position change were large.Between Table 3. List of time periods used for coastal thermo-erosion change detection, bracketed by data acquisitions.Periods referred to as letters correspond to observations covering more than or exactly 1 year.PDD sums were considered as indicators for possible ground ice melt.TD was season-corrected using the number of days with positive mean daily T air occurrence, while TA observations along the cliff bottom did not require correction through open water fraction, because of identical observation period start and end dates.2010 and 2011, the northernmost point of the island retreated by only −4.3 m, while it was −39.4 m during the following period between the two subsequent expedition surveys.However, further away from the cape the opposite picture emerged.Little erosion occurred from 2011 until 2012 when compared to the previous year.Due to these effects, we considered the surrounding 50 m coastline segment as the northern cape area where mean TA over the whole 2010-2013 observation period was −11.6 m a −1 , being still the most rapid erosion along the entire coast of Muostakh Island (Fig. 11, left).Also of note is erosion that occurred between the survey in August 2012 and the GE image, acquired 3 weeks later on 7 September, when the coast next to the northern cape was eroded by up to −7.9 m, because of block failure and collapse due to the deep thermo-niche that existed before in this area (Fig. 14, right), leaving ground ice debris on the sand spit (Fig. 14, left).Despite strong east winds shortly following over 11-13 September 2012 (Fig. 12), TA until July 2013 was only −1.9 m at the cape.

Current thermo-denudation
Using a chronologically consecutive approach of frequent closely spaced TD measurements, we identified seasonal variations for 2010 and interannual variations during the 2010-2013 period.Using intervals of 14-437 days, we dealt with the problem of incomparability of land loss measurements along the cliff top line.The first two consecutive periods 1 and 2 were both 14 days long, while mean absolute TD had doubled from −0.47 ±0.51 in period 1 to −0.98 ±0.52 m in period 2. Period 3 was longer (26 days) and consequently absolute TD was −2.13 ± 0.95 m.Period 4, as the last period which extended over less than a year, absolute TD decreased to −1.84 ± 1.22 m, despite its longer lasting duration of 324 days, 83 days of which have to be considered as the TD active season (50 days in 2010, 33 days in 2011).Period 5 showed mean absolute TD of −5.18 ± 2.88 m over 437 days, 201 days of which nevertheless fell into the TD active season.Period 6 was characterised by very slow TD of −0.3 ± 0.47 m, which is a mixed signal of the remaining 22 days of TD activity in 2012 and 44 days in early 2013.In other words, almost 50 % of the observed absolute TD during the last 3 years happened in 2010.
In order to compare TD intensity across the purely seasonal periods (1-4 and 6), rate calculation of erosion over time through simple annualisation failed, due to large overestimations.For example, simple annualisation would result in mean TD of −29.9 ±13.3 m a −1 during period 3 (summer), while an incredible rapid single value of −79.7 m a −1 would have been measured in period 1 (starting TD in spring) at the northern cape.As a purely descriptive metric, normalisation of TD rates based on specific sums of PDD and OWD for a certain period showed different results for both environmental parameters.OWD-sum-normalised mean TD rates for all periods ranged from −0.2 to −6.6 m a −1 and were on average −2.6 ± 1.4 m a −1 (Fig. 15, right).OWD sum normalisation of TD worked for periods that cover at least one The Cryosphere, 9, 151-178, 2015 sea ice cycle and for period 3, as the only period when TD and TA proceeded unrestricted simultaneously.PDD-sumnormalised rates ranged from −1.7 to −6.3 m a −1 and were on average −3.5 ± 1.6 m a −1 (Fig. 15, middle).The very low standard deviation of σ = 1.6 m a −1 compared to initially σ = 11.5 m a −1 for simply annualised rates, indicates a high degree of levelling between periods and demonstrated TD's dependency on T air .
Based on this finding, we carefully applied a season correction factor into TD rate calculation that accounts for the fraction of number of days with positive mean daily tempera-ture occurrence during the observation period compared with the total annual number (for periods shorter than 1 year) or the annual mean total number (for periods longer than 1 year) of these days (Table 3).Accordingly, the season factor must be ≤ 1.The resulting season-factor-corrected rates were considered to be the actual TD velocity that had taken place.This approach was validated by the fact that the reference annual erosion cycle 2010 (29 June 2010-28 June 2011) turned out to be not affected by the correction approach, since simply annualised and season-factor-corrected TD rates over 2010 were equal, because of a season factor of 1 ( spatio-temporal comparison of TD rates for all periods is shown on maps in Fig. 16. Based on comparable season-factor-corrected rates (Fig. 15, left), we found TD rates during summer were always more rapid (from −8.7 ± 4.6 to −10.2 ± 4.5 m a −1 for periods 2 and 3, respectively) than during spring with −4.1 ± 4.5 m a −1 (period 1) or during fall with −1.4 ± 0.9 m a −1 (period 4) or even −0.2± 0.2 m a −1 (period 6) (Fig. 16).Analyses of interannual variations revealed more rapid TD rates over 2010 of −4.8 ± 2.3 m a −1 , compared to −3.4 ± 1.9 m a −1 during 2011-2012.The 3 years of TD observations, from 2010 to 2013, revealed a mean rate of −3.1 ± 1.6 m a −1 (included in Fig. 15 as reference baseline for current TD), which is a little slower than TA over the same period (Fig. 11,right).Although rapid, it is evident that even on a short-term scale the longer the observation period, the more consistent TD rates will be.
The quality of current thermo-erosion was assessed using the NDTI according to Günther et al. (2012).It turned out that 55 % of the segments from 2010 to 2013 eroded under prevailing TD with an average NDTI of 0.39, while the largest positive TD-indicating NDTI values were observed in the middle of the island, parallel along both the western and eastern coasts with NDTI > 0.5.Along TA-dominated segments, NDTI was −0.34, while on the north-eastern coast, where erosion is most rapid, NDTI was closer to zero, meaning constant erosion of large volumes (Fig. 11,right).
We related positive mean daily T air to TD and TA rates of all observation periods.As a result, we found that a continuous increase of mean daily T air by 1 • C throughout the TD active period is responsible for an acceleration of coastal erosion by −1.2 ± 0.55 m a −1 (Fig. 17).

Open water days
Open water fraction (% d −1 ) was calculated based on SSM/I sea ice concentration data for the past 22 years from 1992 to 2013 in order to identify OWDs and to understand basic characteristics and interannual variability of sea ice in the Buor Khaya Gulf.Background noise present in the data was quantified as 11 ±5.5 % d −1 , by evaluating the open water fraction from December until April, when the land-fast ice zone can be assumed to be completely covered with sea ice, and the sea ice concentration of August, when the coastal waters can be assumed to be free of sea ice.The influence of the perennial Laptev Sea polynya (Reimnitz et al., 1994;Dmitrenko et al., 2005) on the open water data was excluded by taking a 100 km reference zone (Fig. 6).Meier and Stroeve (2008) found SSM/I sea ice concentration data of 15 % matches with the sea ice edge location for the Laptev and East Siberian seas.Adding 11 % d −1 uncertainty due to underestimation, we assumed days with < 26 % within the 100 km to be OWD.
As first-order approximation, we assumed the mean annual open water fraction as indicator for seasonal duration available for TA.Over the past 2 decades, mean open water fraction was 32.1 % d −1 , corresponding to 117 OWD a −1 in the Buor Khaya Gulf, including break-up and freeze-up transition periods.Based on the < 26 % d −1 threshold, the core open water season is on average 81 ± 15 d a −1 long and lasts from 21 July until 8 October (Fig. 18).The difference between 117 and 81 OWD probably reflects 36 days of sea ice drift.The open water season for the current investigation period 2010-2012 was 96 ± 2.5 days long and lasted from 11 July to 14 October.For example Fig. 6 shows Muostakh and the surrounding coastal waters during the very late breakup in 2013.Although the open water season was short in 2013 (73 days), it featured the latest freeze-up (Fig. 18) and was another unusual year, also in terms of mismatch between the TD and TA active seasons (Fig. 12).

Positive degree days
Time series of T air in Tiksi were used to calculate annual PDD sums in Kelvin days (Kd) over the historical time period from 1951 to 2013.The mean seasonal duration based on the first and last occurrence of positive mean daily T air is 133 days, where the season accordingly starts around 17 May and ends on 27 September.However, during this period, days with negative mean daily T air still occur.Mean annual PDD sum over the last 62 a was 660 Kd, the number of days with positive mean T air occurrence was 112 and, accordingly, positive mean daily T air was 5.9 • C. Accounting for the difference of 21 days and assuming a spring-fall partition coefficient of 1 : 3 according to Fig. 19, the season available for TD starts roughly on 1 June and ends on 21 September.Current seasonal and interannual variations of PDD sums from 2010 to 2012 showed that the mean PDD sum was 922 Kd, with a mean daily T air of 7.3 • C, corresponding to seasonal duration of 126 days (Fig. 19).In 2012, the annual PDD sum reached 1010 Kd and exceeded 1000 Kd for the first time in the period of record in Tiksi.A total of 134 days with positive mean daily T air also lengthened the overlap period between the TD and TA active seasons (Fig. 12).

Wind
Wind speed data was cross-checked with data on sea ice concentration and T air for the 2010-2012 study period and for the time series overlap since 1992.During our current study period strong breezes with wind speeds ≥ 10 m s −1 were observed twice during the open water period of 2010 and once in 2011, 2012, and 2013 each (Fig. 12).Severe storm events with wind speeds ≥ 24.5 m s −1 , measuring 10 or higher on the Beaufort wind force scale, usually occurred only during the winter.Over the previous decades (since 1951), during a generalised open water season from 15 July to 15 October, severe storm events in the Buor Khaya Gulf could be .TD during current observation periods from 2010 to 2013 using very high resolution remote sensing data.Symbol size is equivalent to season-factor-corrected TD rates in metres per year (periods 1-6), colour coding expresses proportion between observed seasonal TD rates and mean TD rates over the July 2010-July 2013 period.Seasonal observations show differences of varying TD intensities during short summer (period 3) vs. fall (periods 4 and 6), while interannual observations revealed homogeneous erosion (period 5).
expected to recur at most every 5 years.They had almost exclusively SW direction (Fig. 20), while only < 10 % had northern direction, causing the water level to rise in the Buor Khaya Gulf.Our current study period falls into this storm gap.Generally 6 h wind directions during the TD-active season were diverse, but with prevailing winds from north and north-east and an intensification of strong east winds during the current observation period, which probably enhanced wave action and turbulent heat flux on Muostakh's erosive east facing coast (Fig. 20).Also of note is the increase of T air during the 2010-2012 period that is associated with southerly winds.Despite dominant winds from the north during the sea ice break-up period, usually larger open water fractions and sea ice export out of the Buor Khaya Gulf were associated with winds from the south.In contrast, during the current observation period winds from the north-east and east seem to have favoured sea ice break-up (Fig. 20).

Ground ice and sediment budget
Average macro ground ice content in the subsurface was calculated as 44 ± 4.6 % (n = 1264) by volume.The sediment occupied 56 % of the ground volume.According to Schirrmeister et al. (2003Schirrmeister et al. ( , 2011)), the gravimetric ground ice content of Ice Complex sediments on Muostakh between 0 and 10 m a.s.l. is 108 wt %.Taking this value for the entire Ice Complex series and assuming different densities for ice and solid material according to Strauss et al. (2012), 76 % of the sediment volume is occupied by intrasedimentary ground ice.Combining these 43 % intrasedimentary and 44 % macro ground ice, altogether, 87 % by volume of the geological subsurface is composed of ground ice.
We applied the concept of critical ice content of Are (2012), using the newly available information on topography and ground ice contents.The relative subsidence potential was calculated as δz = 0.69.According to field observations, the ice-poor surface layer was on average 0.7 m thick, the upper two-thirds of which is the active layer.The surface layer was assumed to be mainly composed of organic soil material and therefore has been deducted from the depth of the Ice Complex section.For the northern part of the island, where the base of Ice Complex deposits was detected at 10 m below sea level (Kunitsky, 1989;Grigoriev, 1993) and cliff height is 21 m a.s.l., potential subsidence of the entire 30.3 m Ice Complex thickness equals 20.9 m.If the ground ice would completely melt out, the top of the remaining thawed mate- rial would be situated exactly at sea level.As a possible scenario for the southern part of the island, where no information on the lower Ice Complex boundary is available, we applied the subsidence factor (δz = 0.69) to the mean elevation of the island of 14.4 m, deducted the surface layer, and assumed the Ice Complex base at sea level.Accordingly, potential subsidence is 9.5 m, resulting in a top of thawed material 4.3 m a.s.l.Both scenarios demonstrate that thawing results in a much reduced volume of sediment to be removed from the cliff bottom by waves.In the first scenario, where subsidence would extend almost down to sea level, Are (2012) emphasises that this would mean coastal thermo-erosion may proceed almost exclusively based on its thermal component.
In addition to subsidence, ground ice can result in thermoerosional niches at the water level that undercut the ice-rich cliff.Their theoretical maximum depth is determined by the ice-wedge polygon size (Hoque and Pollard, 2009).Further investigation of the Voronoi diagram was done with a subsample of ice-wedge polygons that are entirely surrounded by other polygons.As a result, each polygon has 5.6 neighbours, on average, which corresponds to the common hexag- Left: wind speed and direction during days with mean daily positive air temperature (T air ), defined as season available for TD to proceed.Right: wind speed and direction during sea ice break-up, defined as time period when open water fraction increases from 30 to 90 %.Middle left and right: mean T air and mean open water fraction depending on prevailing wind direction during TD-active season and sea ice break-up, respectively.Note, although directed away from the erosive coastline, winds from the south also provide favourable conditions for coastal thermo-erosion through higher T air and enhanced sea ice export.onal form of ice-wedge polygons (Christiansen, 2005).They have a mean edge length of 9 ± 4.8 m, occupy a mean area of 162 m 2 , and have a mean diameter of 14.2 ± 2.2 m.Icewedge polygon rows along the cliff bottom and top showed slightly higher macro ground ice contents of 43.5 ± 8.7 and 47.9 ±9.9 vol.%, respectively, compared to 41.1±8.5 vol.% of the middle row (Fig. 21).

Permafrost thaw subsidence
The volumetric content of ground ice in unconsolidated permafrost deposits in East Siberian coastal lowlands exceeds their pore volume in the thawed state (Yershov, 2004), and consequently subsidence results when thawing occurs.The mean annual permafrost thaw subsidence rate of −5.8 ± 2.9 cm a −1 is not randomly or uniformly distributed across the island but varies with geomorphology.In particular, proximity to erosive coastline segments seem to favour thaw subsidence in the island's interior, where for example rates of ≥ near Collinson Head on Herschel Island, a location of constant rapid coastal erosion rates (Lantuit and Pollard, 2008).Remote sensing data showed that seasonal thaw settlement in Alaska is in the range of 1-4 cm (Little et al., 2003) and can be up to 12 cm in places (Liu et al., 2014), while in situ long-term observations of permafrost thaw subsidence are 0.8-1.7 cm a −1 (Shiklomanov et al., 2013).Generally, longterm subsidence occurs by thaw at the top of ice-rich permafrost and drainage of meltwater.Yedoma uplands tend to water drainage and deeper thaw compared to alas depressions (Fyodorov-Davydov et al., 2004).Coastal erosion on Muostakh Island maintains a steep hydrological gradient so that rapid drainage is favoured.Increases in PDD sums result in deeper thaw (Burn, 1998) and permafrost warming (Jonsell et al., 2013).In 2012, we saw the highest PDD sum (1010 Kd) for the entire T air record.Accordingly, ALT on Muostakh in 2012 was on average 47 ± 19 cm.This agrees with a mean thaw depth of 49 cm at the end of August on Samoylov Island in the Lena Delta (Boike et al., 2013).For the nearby Bykovsky Peninsula, where ALT on yedoma uplands varies by climate zone, Fyodorov-Davydov et al. (2008) document a mean seasonal thaw depth of 32 cm during [2003][2004][2005][2006].During these 4 years, the mean annual PDD sum was 670 Kd, very close to the long-term average of 660 Kd but clearly below the average annual PDD sum for the last decade of 775 and 922 Kd for the last 3 years, when ALT on yedoma surfaces increased by 15 cm.
Our comprehensive observations of intensified thaw subsidence in places of shallow ALT (Fig. 10) indicate the opposite to modelling estimations of ALT as a function of surface subsidence by Liu et al. (2012), but are in accordance with the results of ground-based ALT measurements corrected for thaw subsidence by Shiklomanov et al. (2013), and suggest near-surface occurrence of ground ice in areas of strong subsidence.Where deeper thaw encounters ice-rich basal soil horizons or ice-wedge tops, this results not only in active layer deepening, but also in subsidence.Since the upper Ice Complex is ice-rich, increased heat flow into the ground will cause the island to subside.We therefore suggest that the widespread occurrence of peat mounds particularly in the northern part of the island, which present mainly polygonal centres (Fig. 22), is associated with subsidence of the surrounding terrain, rather than with frost heave.
Comparison of our on-site survey with elevation indications from the literature supports the observation of decreasing height of the island.Generally, backshore elevation along the east coast ranges from 21 m a.s.l. in the north to 13 m a.s.l. in the south.Around the former polar station the mean elevation is 6 m a.s.l.Ivanov and Katasonova (1978) report that the height of the island gradually decreases toward the south from 26 to 6 m.Slagoda (2004) presents profiles of Ice Complex sequences that were sampled in 1982 and reached up to 25 m a.s.l.Topographic maps showed the highest point of the island was 25 m a.s.l. in 1982.Grigoriev (2008) presents data of continuous on-site visits since the 1990s and takes 22.6 m as reference height for the north cape of Muostakh.Although all former elevation indications refer to the north cape, it is unlikely that only this portion of the island had an elevation ≥ 25 m a.s.l., but subsidence must have occurred mostly during the last 30 years, particularly in view of the fact that positive deviations from mean PDD sums in NE Eurasia occurred since 1988 (Fedorov et al., 2014).However, our 3-D data set spans 62 years, and estimates of annual subsidence were around −5.8 ± 2.9 cm a −1 over the entire area affected by negative terrain height changes.This is of the same order of magnitude observed elsewhere for ice-rich permafrost (Overduin and Kane, 2006;Fedorov et al., 2012).

Historical and current erosion development
Recently, TA proceeded at −3.4 ±2.7 m a −1 and was therefore 1.9 times more rapid over the past 3 years than the historical record with mean TA of −1.8 ± 1.3 m a −1 .This proportion is consistent with observations made by Günther et al. (2013a), who find recent rates are at least 1.6 times more rapid than historical TD and TA along Ice Complex coastlines throughout the Laptev Sea.They report that for a subset of sites TA accelerated from long-term −3.3 ± 1.2 m a −1 to −5.7 ± 1.2 m a −1 over the past few years.When examined over almost 200 km of Ice Complex coastline and the last 42 years, they find long-term TA was −1.9 ± 1.5 m a −1 , almost identical to Muostakh.This suggests that, despite the annually repeating record-breaking erosion rates on the northern cape, coastal thermo-erosion on Muostakh is not exceptional for Ice Complex coasts in the Laptev Sea.Deviations from the mean can be attributed to local variations of macro ground ice content and exposure to environmental parameters such as to offshore or bay-marine environments.
Since its base area shrank by a quarter (24 %) during 62 years, Muostakh Island will disappear within the next 200 years.If the currently observed erosion rates continue, Muostakh is likely to disappear as island earlier, but no later than within the next 100 years.Examples for disappearing Ice Complex islands on the Laptev Sea shelf exist.According to Gavrilov et al. (2003), the former islands Diomede, Semenovsky, and Vasilievsky have become sandbanks with frozen sediments located very close to the seafloor surface, where thermo-abrasion is still active and proceeds with approximate rates of −0.02 to −0.27 m a −1 .The subsequent submergence of arctic islands results in shoals on the shallow arctic shelf, grounded sea ice pressure ridges (Reimnitz et al., 1994), and loss of island status, as it happened to Dinkum Sands off the Alaskan Beaufort Sea coast (Reimnitz, 2005).Klyuev et al. (1981) report on the shape of Vasilievsky Island in 1823, when the island had a length of 7.4 km and was 463 m wide, quite similar to Muostakh today.As on Muostakh, erosion on Vasilievsky was most intensive along its major axis.Because of its narrow shape, Vasilievsky broke into two parts and was quickly destroyed afterwards due to the unstable ground-ice-rich composition of the island and to the chaotic erosion with rates of up to −100 m a −1 that occurred from all sides until the island was completely destroyed in 1936 (Klyuev et al., 1981).Vasilievsky Island was located quite far offshore and exposed to larger fetch on all sides.Muostakh Island has a narrow central section where we measured TA on both coasts.On the westfacing side, historical and modern rates do not differ greatly (−0.8 ± 0.6 m a −1 ).On the east-facing side, however, modern rates (−2.3 ± 1.5 m a −1 ) are more than twice the historical mean (−1.0 ± 0.2 m a −1 ).Changes to sea ice cover are probably affecting the east-facing coast more, since it is exposed to Buor Khaya Gulf's comparatively large fetch.Possible future break-up of the island will probably occur at this location.Complete disappearance of Muostakh would take away a protection of the Tiksi port harbour approach, but not enhance erosion of the rocky mainland coast.

Interannual and seasonal variability of erosion
Very rapid TD rates of −4.8 ± 2.3 m a −1 in 2010 were followed by slower rates of −3.4 ± 1.9 m a −1 during the 2011-2012 period.Figure 3 suggests that the coastal exposure undergoes geomorphic changes from a high degree of exposed ground ice to an almost complete coverage of ground ice by thermal denudation debris.Remote sensing derived TD is also in agreement with TD from our on-site repeat surveys.According to this control data set TD was −2.7 ± 0.6 m a −1 for 2.8 km of the north-eastern coastline from August 2011 to August 2012.For the entire 2010-2013 cycle we therefore suggest that a constantly high positive mean daily T air in 2010 (Fig. 19) resulted in rapid TD, depositing thawed material and obscuring exposed ground ice in the following year (Fig. 3, top photograph).In 2011, during the very long open water season of 99 days (Fig. 18) TA reworked the coast to predominant steep cliffs (Fig. 23).The main mechanism was block failure due to undercutting by thermo-niches, which probably formed during the strong wind events that occurred in the open water period in 2010 (Fig. 12).However, according to Wobus et al. (2011), thermo-niches develop even under quiet sea conditions.If they have formed during storms, the effect on erosion is combined with a certain delay and consequently not immediately measurable with remote sensing techniques.Maximum thermo-niche development is constrained by ice-wedge polygon size.The continuous presence of thermo-niches was observed in the field exclusively along the north-east and east-facing coast, where backshore height is ≤ 21 m and syngenetic ice-wedges generally stretch across the whole vertical section, indicating that block failure occurs along the vertical plane according to Hoque and Pollard (2009).Given a mean ice-wedge polygon size of 14.2 ± 2.2 m and mean TA of −4.4 ± 2.3 m a −1 in that area, erosion of one complete polygon block would last for 3 years, and our observations captured this time frame within the 2010-2012 period.Around the northern cape, current TA was −11.6 m a −1 (the most rapid segment in Fig. 11), so that erosion of one complete polygon block takes roughly 1 year there.In general, we suggest that coastal erosion rates are linked to geomorphology particularly through the varying macro ground ice content and its subsurface distribution.
The established dependency of coastal thermo-erosion on T air (Fig. 17 +4 • C, which also seems to be a threshold for TD activity on Muostakh (Fig. 17).In West Siberia, where massive ground ice occurs along the arctic coast, Vasiliev et al. (2006) find that coastal retreat was 2 times higher where ground ice content was 45 %, when compared to places with only 25 %.Lantz and Kokelj (2008) studied retrogressive thaw slumps in the Mackenzie Delta region and find that an increase in mean summer T air by 1.3 • C leads to a 1.4 times more rapid general slump growth and an almost doubling of the slump headwall retreat.Our result of currently 1.9 times more rapid erosion is also consistent with recently observed 1.6 times more rapid coastal erosion rates in the Laptev Sea region by Günther et al. (2013a).Wobus et al. (2011) used time-lapse photography to study thermo-erosion along Alaska's Beaufort Sea coast.They report thaw rates of 1-6 cm d −1 during spring, prior to sea ice break-up.During our periods 1 and 2 (15 June-13 July 2010), both lasting for 14 days (Fig. 16, top left and top middle), absolute cliff top line retreat was 0.47 and 0.98 m, and mean daily T air was 6.1 • C and 10.6 • C, respectively (Table 3).With respect to largely exposed ground ice, this corresponds to thaw rates of 3-7 cm d −1 or ablation rates of 5.5-6.6 mm d −1 • C −1 , according to Braithwaite (1995).Wobus et al. (2011) 2012) turned out to be sensitive primarily to meteorological parameters, not to sea ice extent.Although we did not account for the important energy supply through solar irradiation during sunny and overcast days (Lewkowicz, 1986), nor for sea surface temperature, our TD normalisation efforts also showed TD sensitivity to T air rather than to OWD (Fig. 15).Future studies should extend seasonal analyses to TA, to better account for marine forcing.

Changes in environmental parameters
The recent lengthening of the sea-ice-free season by 15 days from 81 (1992-2013 mean) to 96 days on average for the 2010-2012 erosion observation period was due to 10 days of prolongation in early summer and 5 days of later freezeup within the Buor Khaya Gulf.This is in accordance with a general trend of accelerating early ice retreat in the Laptev Sea, probably because of thinner sea ice (Krumpen et al., 2013).Based on our observations, on average, the open water season starts on 20 July.Despite the fact that interannual variations are considerable (35 OWD in 1996vs. 99 OWD in 2011), Karklin and Karelin (2009) report that negative anomalies of earlier break-up became a rule from 1999 to 2005.Our data show a continuation and strengthening of this trend (Fig. 18), parallel to rising T air (Fig. 13).This also adds perspective to the trend of OWD increase of 0.5-1 d a −1 observed by Barnhart et al. (2014) for coastal cells in the Laptev Sea with significant trends over the entire 1979-2012 SSM/I data set.Overeem et al. (2011) show that open water duration has been increasing along the Alaskan Beaufort Sea coastline from icy to ice-free over decades based on SSM/I ice cover calculations.In our case, most likely, the warm Lena River waters are likely to additionally support local seasonality via early break-up, because of generally thinner sea ice (Spreen et al., 2011) and earlier spring floods (Federova et al., 2009).Even around Muostakh, spatial variations of sea ice conditions are pronounced, where the rapidly eroding east-facing coast can be completely free of ice, while on the west-facing coast, towards the protected Tiksi Bay, sea ice may persist.
According to Gukov (2001), the waters west of the island are ice-covered for 12 days longer than in Buor Khaya Gulf.However, depending on either thermal or mechanical sea ice break-up (Petrich et al., 2012), the opposite picture might emerge, as it was in 2013 (Fig. 6).
Wind speeds at Tiksi range between 0 and 25 m s −1 .Strong winds exceeding 10 m s −1 occur almost exclusively during the winter months, with continuous sea ice cover (Fig. 12).Adding to this seasonal protection of the island from wind-driven wave action is the fact that almost all high wind speed events come from the south-west to south direction (Fig. 20), causing water level to fall.Even when ice-free, the maximum fetch south to south-westward of Muostakh is less than 50 km.In terms of constant heavy swells, the earlier start of the open water season might have strong impacts on TA on Muostakh Island, since winds during the sea ice breakup period are almost exclusively from the north and northwestern directions, where Muostakh is exposed to nearly unlimited fetch, causing water level to rise.
The very warm summer of 2010 with a positive mean daily T air of 7.7 ±5.3 • C, during 124 days with positive mean T air occurrence, was also accompanied by high net radiation (Boike et al., 2013).In 2011, T air was lower (6.5 ± 4.9 • C over 121 days) but winds were stronger, enhancing the convective heat transfer and maintaining melt on ground ice exposures.According to Langer et al. (2011), wind speed in the Lena Delta features a diurnal pattern, with enhanced heat exchange during the day and lowered turbulence during the night.In the very warm summer of 2012, positive mean daily T air was 7.5 ± 4.7 • C over the 134-day long TD-active period.Mean July T air of 12.5 • C in 2010, 10.7 • C in 2011, and 11.5 • C in 2012 were higher than the long-term mean July T air of 7.5 • C (1951-2012).Thus, not even the criterion for a subpolar climate, that T air is below 10 • C in the warmest month (Neef, 1956), was met over these 3 years.Besides this warming in the tundra zone, Fedorov et al. (2014) report that the greatest recent summer warming was observed in the adjacent forest tundra and northern taiga.

Macro ground ice as local controlling factor
Our estimates of 87 % total volumetric ground ice content on Muostakh are slightly above previous specifications of 80 % (Are, 1988a;Slagoda, 2004).Our polygon mapping approach allows sediment centres to touch the polygon's boundary and accordingly mapped ice wedges may have zero width in places.Thus, our geomorphometric method probably overestimates polygon sediment centre size, resulting in conservative rather than excessive assessments of macro ground ice content.Since the sediment on Muostakh is one of the most coarsely grained and poorly sorted examples of Ice Complex (Slagoda, 2004;Schirrmeister et al., 2011), baydzharakhs there are better preserved following thawing, in contrast to those elsewhere which quickly slump or undergo transport (Are et al., 2005).
The shape of syngenetic ice wedges in a vertical plan often deviate from the ideal wedge form (Popov et al., 1985).Generally, 2-4 rows of baydzharakhs could be observed on the slope.Our results imply that macro ground ice distribution has a non-uniform vertical hourglass shape, with higher ice contents and ice wedge widths at the top and bottom of the slope, and probably below sea level.Variable sedimentation and preservation conditions during the Late Pleistocene defined Ice Complex accumulation rates (Wetterich et al., 2011).Thaw unconformities and several ice-wedge generations that are nested in one another also caused varying ice-wedge width across the profile.Field observations also confirm that Pleistocene syngenetic ice wedges on Muostakh Island are thickest at the bottom (rather constant width and not increasing width) and at the top of the coastal exposure, where at the upper 4 m of the profile Holocene ice wedges are often interposed.The higher ice contents thus occur at the positions where we measure coastal erosion and are predisposed to favour intense coastline retreat as a result of warming (TD) and wave action (TA).

Conclusions
In this study, we found that continuous coastal erosion on Muostakh Island in the Buor Khaya Gulf of the Laptev Sea during the last 62 years caused the land loss of about 24 % of the island's area, while its volume shrank by 40 %.Muostakh is composed of Ice Complex permafrost deposits, of which up to 87 % of the subsurface is occupied by ground ice.This exposes the island to thermal disturbances from coastal erosion and seasonal thaw and permafrost thaw subsidence, leading to further destruction of the island and its final disappearance, which we expect to take place within the next 100 years under recently changing environmental conditions.Average subsidence of the island's surface was −3.6 ±1.8 m over 62 years and therefore −5.8 ±2.9 cm a −1 .Average coastal erosion from 1951 to 2012 was −1.8 ± 1.3 m a −1 , while recent rates were 1.9 times more rapid at −3.4 ± 2.7 m a −1 .At the highly erosive northern cape of the island, the distance from the mainland has been increasing by −9.6 m a −1 ; this value also accelerated during the recent past to −17 m a −1 .
Our systematic seasonal analyses demonstrate that the currently higher intensities of the two coastal erosion processes, thermo-abrasion (TA) and thermo-denudation (TD), are controlled at least in part by the increasing open water season and summer air temperatures (T air ).The open water season has lengthened from 81 open water days (OWD) on average for the past 2 decades by 15 OWD over the 2010-2012 observation period and, for example, up to 99 OWD in 2011.Annual positive degree day (PDD) sums from 1951 to 2013 in the nearby town of Tiksi were 660 Kd on average and strongly increased to 1010 Kd in 2012.Accordingly, the seasonal duration available for TD has lengthened from 110 days on average to 134 days in 2012.
We show that normalisation of diverse TD rates (σ = 11.5 m a −1 ) through PDD sums for each period decreases variability of TD rates across all subperiods to σ = 1.6 m a −1 .We therefore propose that TD rates for short season time periods are corrected by a seasonal factor to make them comparable.Our observations suggest that coastal erosion on this Ice Complex coast increases by 1.2 m a −1 per 1 • C of positive mean daily T air .
Interannual variations of coastal erosion are also related to local factors, such as macro ground ice content variation with depth, affecting thermo-niche development, which depends on ice-wedge polygon size.
We found a phase shift between TA and TD.TA is only active during the open water season, while TD can proceed throughout the summer.In June, when T air rises and the mainland is already free of snow, mud flows from TD accumulate on top of snow at the cliff bottom, while sea ice is preventing TA.In late August, TD is slowed due to refreezing of the active layer and coastal slopes, while TA may still proceed until land-fast sea ice develops.Currently observed shifts in sea ice duration and T air result primarily in an extension in overlap time of the active seasons for TA and TD -the resulting simultaneity of both processes is more important than the extension of either active season, because the widely varying value domains of TD and TA align and mutually reinforce each other to more effective erosion of thawed permafrost and associated mass displacement.

Figure 1 .
Figure 1.Left: situation of the right-hand map in East Siberia (Russia; source: ESRI).Right: location of Muostakh Island within Buor Khaya Gulf, central Laptev Sea (September 2010 Landsat-5 imagery as background)

Figure 2 .Figure 3 .
Figure 2. Coastal thermo-erosion over time in the northern part of Muostakh Island.Historical cliff bottom (1951, blue) and current cliff top line (2012, red) border the subaerial coastal thermo-erosion zone.Left: orthophoto of historical aerial imagery draped over 1951 DEM.Note: the lake in the lower left has been drained during the observation period.Right: oblique photograph taken from helicopter in August 2012.Location of 2011 and 2012 field camp, maximum island width for scale.

Figure 4 .
Figure 4. Map showing locations of active layer thickness measurements on Muostakh made during field work in August 2012.Sampling points were classified and colour coded in steps of 10 cm.2012 GeoEye image as background.

Figure 5 .
Figure 5. Elevation dependent bias in 2013-1951 raw DEM difference and correction using an empirical function of second order polynomial.

Figure 6 .
Figure 6.Detail of a Landsat-8 image (5, 4, 2 CIR) showing sea ice break-up around Muostakh Island (red square, size corresponds to one 12.5 km SSM/I Pixel) on 18 July 2013.The open water fraction on this day, 100 km around Muostakh (yellow circle), was 20 %.The sea-ice-free season started 19 days later on 6 August 2013.

Figure 7 .
Figure 7. Top: photograph of a baydzharakh field on the east coast of Muostakh Island.Bottom: examples of slope profiles across baydzharakhs, showing differences in baydzharakh spacing.

Figure 8 .
Figure 8. Left: points mark mapped baydzharakh centre locations, used for derivation of the Voronoi diagram.Middle: determination of the Euclidean distance within each polygon.Right: construction of largest possible circles within polygons using maximum Euclidean distance as radius, representing the sediment component of the subsurface.Calculation and interpolation of macro ground ice content between circles based on the ratio of area occupied by circles and total polygon area (8 August 2010 WorldView-1 imagery as background).

Figure 9 .
Figure 9. Left: DEM from 1951 stereoscopic aerial photography.Middle: DEM from 2013 GeoEye stereo pair.Symbol size is the classified planimetric coastal erosion rate.Colour-code displays volumetric erosion from 1951 until 2010 for 118 coastline segments.Right: difference raster from multitemporal DEMs representing elevation changes over 62 years.

Figure 10 .
Figure 10.Classified observations of active layer thickness (ALT) in August 2012 in relation to mean annual permafrost thaw subsidence from 1951 to 2013 show intensified subsidence in places of shallow ALT.The outliers are indicated by open circles.The colour code is adjusted to ALT classification in Fig. 4. Inverse relation of ALT and subsidence indicates water drainage on the permafrostactive layer interface and consequently irreversible ground ice thaw.
The 13 July 2010 image represents the T air summer peak and marks the start of the open water season.The fourth image was acquired during the open water season, when mean daily T air had already begun to fall (8 August 2010).Together with the previous image, it completely spans a period when both TD and TA are active.The 2010 fall season was bracketed by the fourth and the fifth image, acquired in early 2011 during ice melt and rising T air .The 7 September 2012 image, acquired at the peak of the open water season and falling T air , together with the previous image, captured almost two complete seasonal cycles.The 17 July 2013 image was acquired prior to late sea ice break-up in 2013, and completes not only fall 2012, but captured also spring 2013.

Figure 11 .Figure 12 .
Figure 11.Point by point coastal thermo-erosion over time, red crosses indicate mean values.Left: historical thermo-abrasion (TA) from 1951 to 2012 in relation to current changes from 2010 to 2013, showing TA acceleration at 95 out of 118 coastal segments.Line of equal rates.Right: current thermo-denudation (TD) in relation to TA, showing TD-and TA-dominated coastal erosion regimes were around the same frequency.Virtual normalised difference thermo-erosion index (NDTI) zero-line for differentiation.

Figure 13 .
Figure13.Open water fraction and positive degree-days (PDDs) as annual sums for the reference period of both records.

Figure 14 .Figure 15 .
Figure 14.Close up view of coastal erosion at the northern cape of Muostakh Island.Left: selected cliff top position lines of the 2010-2013 period.Blue line outlines cliff bottom position in early 2010.Ends of cliff top lines are at sea level, marking the northernmost point of the island.Note stranded ground ice debris on the sand spit.Right: photograph of the northern cape (8 August 2012).
Figure16.TD during current observation periods from 2010 to 2013 using very high resolution remote sensing data.Symbol size is equivalent to season-factor-corrected TD rates in metres per year (periods 1-6), colour coding expresses proportion between observed seasonal TD rates and mean TD rates over the July 2010-July 2013 period.Seasonal observations show differences of varying TD intensities during short summer (period 3) vs. fall (periods 4 and 6), while interannual observations revealed homogeneous erosion (period 5).

Figure 17 .Figure 18 .
Figure 17.Mean annual coastal thermo-erosion (eight data points for TD and three for TA) vs. positive mean daily T air .Each erosion data point is a mean of all 118 coastal sections; error bars indicate standard deviation of rates.T air was calculated for specific observation periods using only days with mean T air >0 • C.

Figure 19 .
Figure 19.Time series of positive mean daily T air over the period from 1951 to 2013 show lengthening and intensification.Seasonal duration available for thermo-denudation lengthened from 111 days for the entire period to 127 days on average for complete years of the current observation period (2010-2012).Simultaneously, positive mean daily T air was currently 1.4 • C higher and increased from 5.9 to 7.3 • C.

Figure 20 .
Figure 20.Wind charts and rose diagrams for the long-term reference time period (top row) and the current observation periods (bottom row).Left: wind speed and direction during days with mean daily positive air temperature (T air ), defined as season available for TD to proceed.Right: wind speed and direction during sea ice break-up, defined as time period when open water fraction increases from 30 to 90 %.Middle left and right: mean T air and mean open water fraction depending on prevailing wind direction during TD-active season and sea ice break-up, respectively.Note, although directed away from the erosive coastline, winds from the south also provide favourable conditions for coastal thermo-erosion through higher T air and enhanced sea ice export.

Figure 21 .
Figure 21.Histograms of macro ground ice content along erosive coastal cliffs on Muostakh Island.Classification of ice-wedge polygons into three different vertical positions shows shift towards higher macro ground ice contents for the lower and upper parts of the subsurface.Vertical lines indicate mean values.
et al.: The disappearing East Siberian Arctic island Muostakh −11 cm a −1 were observed close to the northern cape, eroding at ≤ −10 m a −1 .Similarly, Short et al. (2011) showed −10 to −15 cm of terrain displacement during summer 2010

Figure 22 .
Figure 22. of peat mounds (incipient baydzharakhs) in the northern part of Muostakh, where the island's surface is affected by strong permafrost thaw subsidence (person for scale).

Figure 23 .
Figure 23.Photograph of catastrophic collapse through undercutting, reflecting activation of TA during the rising water level at the end of August 2011, when sea water comes in contact with permafrost at the level of thermo-niches.Under calm weather conditions the ground ice block decayed within 2 days.
) and the seasonally varying intensity of thaw on coastal ice cliffs agrees with other observations.Grigoriev et al. (2006) link coastal dynamics of East Siberian coasts to changes in climatic and permafrost conditions.They find a positive correlation between coastal retreat rate and macro ground ice content, especially when T air exceed www.the-cryosphere.net/9

Table 1 .
List of very high resolution satellite imagery used for change detection and summary of multisensor bundle block adjustment.

Table 2 .
Volumetric losses and associated mass displacement on Muostakh Island, based on DEMs of 1951 and 2013 for different compartments of the subsurface, assuming fractional volumes of 44 % macro ground ice and 43 % intrasedimentary ground ice.

Table 3 ). The www.the-cryosphere.net/9/151/2015/ The Cryosphere, 9, 151-178, 2015 F. Günther et al.: The disappearing East Siberian Arctic island Muostakh
also observe acceleration after open water season begins, which is consistent with more rapid rates on Muostakh during period 3 (13 July-8 August 2010), the only phase when TD and TA proceeded simultaneously, and thaw was on average 8 cm d −1 and, in places, reached 17 cm d −1 under a mean daily T air of 12.7 • C. Thaw was slowest during spring 2013 at 0.4 cm d −1 .The niche and block erosion model ofRavens  et al. (