Coastal Erosion Variability at the Southern Laptev Sea Linked to Winter Sea Ice and the Arctic Oscillation

Arctic coastal erosion experiences pronounced effects from ongoing climate change. The Laptev Sea figures among the Arctic regions with the most severe erosion rates. Here, we use unprecedentedly long records of almost 30 years of annual in‐situ coastal erosion rate measurements from Bykovsky Peninsula and Muostakh Island to separate the main modes of variability, which we attribute to large‐scale drivers. The first (lower‐frequency) and second (higher‐frequency) modes are associated with winter sea‐ice cover in the Laptev Sea and with the Arctic Oscillation, respectively, which together account for 85.1±24.1% of the total observed variance. Arctic coastal erosion has so far been neglected in Earth system models (ESMs). The proposed mechanisms set favorable conditions for coastal erosion at large scales (synoptic to planetary scales), compatible with those represented in modern ESMs.


Introduction
The Arctic has been experiencing pronounced effects of climate change: increasing surface (Serreze et al., 2009) and permafrost temperatures (Biskaborn et al., 2019) and decreasing sea-ice extent (Notz & Stroeve, 2016). Consequently, the Arctic coasts are now being longer than previously exposed to the action of warm air and ocean waves, leading to thermal and mechanical erosional processes. At the Laptev Sea, eastern Siberian Arctic, the historical mean coastal erosion rate is estimated at 0.7 m year −1 , somewhat larger than the pan-Arctic mean of 0.5 m year −1 (Lantuit et al., 2012), with specific locations showing annual retreats of >20 m  and figuring among the most rapidly eroding sites in the Arctic. Irrgang et al. (2018) presented similar mean rates of 0.7 m year −1 with significant time variability for the Yukon Coast, Canada.
The erosion of permafrost coasts release substantial amounts of organic carbon (OC) to the marginal Arctic seas, estimated at 14 ± 4 Tg year −1 , comparable to the pan-Arctic OC flux from riverine discharge of 40 ± 4 Tg year −1 (Wegner et al., 2015). Vonk et al. (2012) estimated the OC flux from the East Siberian Arctic Shelf (comprehending part of the Laptev and East Siberian Seas) at 44 ± 10 Tg year −1 , from which 11 ± 4 Tg year −1 would be resulting from coastal erosion. Günther et al. (2013) estimated a significantly smaller carbon flux of 0.66 ± 0.05 Tg year −1 from the Laptev Sea coast, using a combination of remote sensing data and in-situ measurements. Moreover, the organic matter from thawed permafrost was found to be highly biologically reactive (Vonk et al., 2014). Couture et al. (2018) estimated that only 13% of the eroded OC is sequestered in the nearshore sediment at the Beaufort Sea. In an incubation experiment, Tanski et al. (2019) recently showed that the direct release of atmospheric CO 2 from coastal erosion has been underestimated. Given the projected rapid decrease in Arctic sea-ice extent (Barnhart et al., 2016), intensification of surface wind and waves (Dobrynin et al., 2012(Dobrynin et al., , 2015 and permafrost thaw (Schaefer et al., 2014;Schuur et al., 2009), the degradation of ice-rich thermo-abrasive Arctic coasts will likely increase in the future, leading to the recycling of larger amounts of permafrost carbon into CO 2 and, thus, increasing the atmospheric carbon-climate feedback. However, comprehensive modeling studies considering the role of coastal erosion on climate are still needed (Fritz et al., 2017).
Arctic coastal erosion is primarily limited by the presence of sea ice. Overeem et al. (2011) suggested that the duration of the open-water season (OWS), that is, the time of year when the coast is sea-ice free, is a good first-order estimator for coastal erosion, given that both quantities have increased at similar rates between 1979-2002 and 2002-2007 (1.5-fold and 1.6-fold, respectively).  presented a circum-Arctic analysis of OWS duration and concluded that it has larger relevance in explaining coastal vulnerability, in comparison with setup and wave heights. Both studies used sea-ice concentration (SIC) data from nearshore grid cells to derive OWS duration and disposed of long-term means of coastal erosion rates. The variability of Arctic coastal erosion also responds to environmental conditions within the OWS, such as positive air surface temperatures  and the frequency and intensity of storms (Jones et al., 2009). Cunliffe et al. (2019) recently showed that single-storm events may be relevant for total annual shoreline changes.
The Arctic Oscillation (AO) is the dominant mode of sea-level pressure variability in the Northern Hemisphere (Thompson & Wallace, 1998) and plays an important role on weather, including the frequency and intensity of storms (Thompson & Wallace, 2001) and Arctic SIC anomalies (Wang & Ikeda, 2000). During winter, its positive phase (AO+) favors positive surface air temperature anomalies over the Eurasian continent (Thompson & Wallace, 2000). At the Laptev Sea, the winter AO+ is associated with surface winds blowing from southeast, its continental margin, transporting sea ice away from the coast and into the central Arctic Ocean, thus contributing to the opening of polynyas and formation of new thin ice, increased heat fluxes from ocean to atmosphere, and local surface warming (Rigor et al., 2002). Krumpen et al. (2013) and Itkin and Krumpen (2017) showed that late-winter sea-ice export from the Laptev Sea, thus the formation of thin ice, is correlated with negative SIC anomalies in the forthcoming summer and may as well contribute to earlier onsets of the melt season.
Previous studies have focused on process-based approaches to address the issue of Arctic coastal erosion by, for example, modeling block-failure events (e.g., Hoque & Pollard, 2009;Ravens et al., 2011). Although physically meaningful, the high spatial resolution needed in that setup (order of meters) is not compatible with the scale of the still relatively coarse-resolution state-of-the-art Earth system models (ESMs, order of hundreds of kilometers). Here, we aim at exploring the predominant large-scale drivers of coastal erosion observed at the southern Laptev Sea, by encompassing information from large areas (synoptic to planetary scales), thus responding more directly to dynamic and thermodynamic mechanisms of the climate system, which are inherently better represented in modern ESMs than small-scale ones.

Data and Methods
The search for statistically robust relationships between external drivers and coastal erosion variability is often hampered by the lack of long and well-resolved observations (Lantuit et al., 2011). Although in the last decade, the availability of high-resolution satellite imagery has advanced the use of remotely sensed shoreline-change mapping (Jones et al., 2018). Here, we analyze unprecedentedly long in situ observations of coastal erosion rates from Bykovsky Peninsula and Muostakh Island at its North Cape (hereafter, Muostakh-N) and at its northeast coast (Muostakh-NE), southern Laptev Sea (Figures 1a and 1b). Measurements have been yearly made by the end of August since 1982 until present (Grigoriev, 2019). Some gaps in observations occur between 1983 and 1996 in different years for each site. Since 1982, a total of 28 years of data are available for Bykovsky and Muostakh-NE and 29 years for Muostakh-N. Since the observations are spatially limited, generalizations to the scale of the Laptev Sea must be made with caution (see Location of the long-term key monitoring sites. The background satellite image (9 September 2009) was obtained from NASA Worldview application (https://worldview.earthdata.nasa.gov) and had contrast enhanced. Power spectral density (PSD) of (c) FMA SIC averaged over the Laptev Sea (FMA LSIC) and of (e) FMA and (f) JJA AO indices. The green full and dashed lines indicate the 90% and 95% confidence levels, respectively. In (d), we show a map of predominant periodicity (years), corresponding to the maxima of PSD in FMA SIC per grid cell. Note that the frequencies are obtained from the 39-year-long period  from ERA-Interim. Even though peaks in PSD are significant at the 95% level, calculated with a formal probabilistic method (see the supporting information), the exact periods should be taken with caution due to the frequency discretization. section 4 for a detailed discussion). For more detailed descriptions of the key monitoring sites, we refer the reader to Lantuit et al. (2011) on Bykovsky Peninsula and Günther et al. (2015) and Overduin et al. (2016) on Muostakh Island. The main modes of variability of coastal erosion observations are separated with a principal component analysis.
The AO is defined as the first empirical orthogonal function, of SLP north of 20 • N ( Thompson & Wallace, 1998). The AO index is the principal component, associated with the first empirical orthogonal function. We calculate the AO in seasonal means and focus on the Arctic winter (February-March-April [FMA]), when sea ice reaches its maximum concentration in the Laptev Sea, and in summer (June-July-August [JJA]), before the field measurements.
From ERA-Interim reanalysis (Dee et al., 2011), we obtain data on sea-ice concentration, SIC (%), for the calculation of seasonal means and OWS duration, sea-level pressure, SLP (hPa), for the calculation of the AO indices, and a list of dynamics and thermodynamics variables to explore the underlying physical mechanisms (see the supporting information). These data are disposed in a Gaussian grid of 0.7 • horizontal resolution, from 1979 to 2018. Daily SIC data are averaged over the Laptev Sea (100 Figure 1a) to create Laptev SIC (LSIC) time series. From daily LSIC, we use a threshold of 15% open ice (85% sea-ice cover) to determine the start date of the melting period, the end date of the recovery period, and the duration of the OWS (see the supporting information for details).
We define especially strong and weak erosion rates, those larger than half a standard deviation above or below their mean (>|0.5 |, "extreme"). Analogously, we define close-to-neutral conditions those when coastal erosion rates were smaller than half a standard deviation around their mean (<|0.5 |, "neutral"). In order to focus on the interannual variability, long-term linear trends are removed from all time series, and anomalies were calculated with respect to the 1979-2018 period. More details on the statistical methods and metrics are available in the supporting information.

Laptev Sea Ice and the Lower-Frequency Mode
Winter (FMA) SIC presents predominantly low-frequency variability over much of the eastern marginal Arctic seas. The LSIC varies in the decadal timescale, with a peak period of ∼19 years (Figure 1c), similarly to its neighboring seas. Spatially, most of the Laptev, East Siberian, and Chuckchi Seas FMA SIC shows frequency with maximum power spectral density values corresponding to periods longer than 15 years (Figure 1d), as well as in parts of the East Kara Sea, especially in non-coastal grid cells, suggesting independence from land processes. In the Barents and Beaufort Seas, variability of FMA SIC is relatively higher, with peak periods shorter than 10 and 5 years, respectively.
The predominant low-frequency variability of FMA LSIC is likely due to oceanic mechanisms. Zhang (2015) showed that Atlantic and Pacific heat transport have maximum significant correlations with September Arctic sea-ice extent with a 2-year lag. The Atlantic Multidecadal Variability (AMV), in its positive phase (AMV+), is associated with a negative SLP anomaly over the Kara, Laptev, and East Siberian Seas, driving a cyclonic circulation and consequent export of sea ice into the Arctic Ocean during winter Castruccio et al. (2019). The AMV+ is also significantly correlated with positive anomalies of cloud longwave radiative forcing at the surface, consequent positive surface air temperature anomalies, and negative anomalies of ice strength during winter over much of the East Siberian and Laptev Seas (Castruccio et al., 2019). Although here noted, the study of the low-frequency variability source of winter LSIC is beyond the scope of this letter.
Negative winter LSIC anomalies are associated with longer OWS (r = −0.57 in March, p < 0.05, Figures 2a and 2d) and with larger coastal erosion rates (r < −0.50, p < 0.05, Figure 2g). Moreover, positive FMA LSIC anomalies are more strongly associated with late onset dates of OWS (r = 0.50, p < 0.01), than with early demise dates of previous year's OWS (r = −0.33, p ≃ 0.21). Indeed, late-winter sea-ice export from the Laptev Sea persists into negative summer SIC anomalies (Krumpen et al., 2013), which anticipates the onset of the melting period (Itkin & Krumpen, 2017). Thus, negative FMA LSIC anomalies contribute to lengthen the OWS, which has been suggested to be the first-order driver of Arctic coastal erosion Overeem et al., 2011).
The first principal component (PC1), derived from coastal erosion observations, accounts for 56.5 ± 16.0% of the total variance and seems to be associated with FMA LSIC anomalies (Figure 3a). PC1 shows minimum  Figure S1). Applying a low-pass filter to PC1 as well, a linear relationship with FMA LSIC and the OWS duration appears even stronger (up to r = −0.89, Figure S2). Therefore, we propose that PC1 represents the low-frequency variability in FMA LSIC, which modulates the duration of the OWS.

The AO and the Higher-Frequency Mode
In parallel, coastal erosion seems to be associated with the AO in both winter and summer, with opposite polarities. The differences between the strong and weak erosion composites of SLP show spatial patterns strikingly similar to the AO in its positive phase (AO+) in FMA (Figures 2b and 2e) and in its negative phase (AO−) in JJA (Figures 2c and 2f). Not only is the general annular mode visible in the composite means of SLP and Z200, but also is the meridional shift of the node latitude between winter and summer (Ogi et al., 2004), suggesting an association between coastal erosion and the AO with its equivalent barotropic signature, in accordance with Thompson and Wallace (1998).
The second principal component (PC2) from coastal erosion observations accounts for 28.6 ± 8.1% of the total variance and varies in quasi-unison with the AO indices, especially after 2010 (Figure 3c). Correlation coefficients are r = 0.46 (p < 0.05) with the FMA AO and r = −0.33 (p = 0.25) with the JJA AO. Bykovsky and Muostakh-N present large PC2 loadings with opposite signs, in accordance with the anti-correlation found between the FMA and JJA AO in years of especially strong and weak erosion rates, explored in the next section. To support the PC2-AO linkage, we propose the following mechanisms.
The FMA AO is associated with southern surface winds, which advect sea ice away from the coast of the Laptev and East Siberian Seas, increasing the production of new thin ice, the opening of polynyas and flaw leads, thus increasing heat fluxes from the ocean to the atmosphere, and resulting in local surface warming (Rigor et al., 2002). We also identify this pattern in regressions using the FMA AO index and in our PC2 erosion composite difference (Figures S3a and S3b).
The JJA AO− is associated with positive surface temperature anomalies in the Laptev Sea. Ding et al. (2017) showed that anticyclonic circulation anomalies with barotropic structure over Greenland the Arctic Ocean (the JJA AO polar center of action, Figure 2c) increase low-level cloudiness by warming and moistening the lower troposphere, which is then translated in increasing incoming longwave radiation. They further suggest that this mechanism is dominated by circulation changes, reinforcing the role of the AO. We verify this mechanism in PC2 composite differences ( Figures S4c-S4h). We find significant correlations between the JJA AO index and surface downward longwave radiation over the Laptev Sea. Indeed, significant regression coefficients seem bounded by the coastline, suggesting that the Laptev Sea plays an important role as a source of moisture to this end. However, total integrated water-column content does not show such significant relationships, despite of a pattern of predominant positive anomalies. The signal is made more clear in low cloud-cover anomalies, suggesting that moisture is mainly increased in the lower troposphere. Therefore, a negative JJA AO condition would, in a large scale, create favorable conditions for thermally driven coastal erosion by increasing surface warming. The JJA AO− is also associated with negative summertime Arctic SIC anomalies, driven by negative SLP anomalies over the polar cap and an associated anticyclonic circulation (Ogi et al., 2008;Wernli & Papritz, 2018). In fact, Ogi et al. (2016) showed that both winter and summer AO indices are good predictors for September SIC anomalies. Negative SIC during summer, and consequent larger fetch for wind waves, would also increase the vulnerability of Arctic coasts to the action of waves.
While PC1 and FMA LSIC show variability in decadal timescales, the AO shows relatively higher-frequency variability with predominant periods between ∼2.7 and ∼4.2 years, for both winter and summer indices (Figures 1e and 1f). Therefore, we suggest that PC2 represents the higher-frequency imprint of the AO onto the coastal erosion observations. We do not attribute the third principal component to any climatic driver; it probably represents the signature of local properties, not responding directly to the FMA LSIC and AO effects (Figure 3e). The first two modes account together for 85.1 ± 24.1% of the total variance and can be associated with mechanisms (see summary in Figure S6) identifiable at scales larger than those typically considered for Arctic coastal erosion.

Quantifying the Role of the Large-Scale Drivers at Each Site
We employ multiple linear regression (MLR) models to estimate the relative role of FMA LSIC and the AO indices at explaining the variance of observed coastal erosion rates at each individual site. To emphasize the different frequencies between drivers, a low-pass Lanczos filter with a cutoff-period of 10 years is applied to FMA LSIC, allowing decadal to longer-scale variability.
The role of each explanatory variable varies considerably between sites. The FMA LSIC alone explains 21%, 30%, and 22% of the observed variance of coastal erosion at Bykovsky, Muostakh-N, and Muostakh-NE, respectively. The AO indices explain at most 17% (FMA) and 27% (JJA) of the variance at Bykovsky. The proportion of explained variance increases, but does not double, when both AO indices are used as covariates, suggesting some degree of collinearity. Taking the full model, that is, combining FMA LSIC and both FMA and JJA AO indices, 53% (Bykovsky), 28% (Muostakh-N), and 20% (Muostakh-NE) of the variance is explained. MLR experiments are also done separately taking only years of extreme and neutral erosion rates. The goodness-of-fit of models increases in 77% of the experiments in the extreme case and decreases in 72% of the experiments in the neutral case. Thus, the proposed drivers are generally better able to explain coastal erosion in years of extreme values, than during years when coastal erosion rates are closer to their mean state. Therefore, the MLR models capture not only the low, but also the high values of erosion rates, suggesting that the proposed drivers are associated with the anomalies at both ends of the observed range. For Bykovsky, the proportion of explained variance reaches 87% for the full model with extreme years (r = 0.95). The single-variable FMA LSIC model performs best in the extreme cases and explains 53% of the variance for Muostakh-N (r = 0.75) and 32% Muostakh-NE (r = 0.60). The dominant role of the low-frequency FMA LSIC variability is visible in the three modeled time series (Figures 3b, 3d, and 3f).
The differences in magnitude in MLR coefficients (Table S1) Figure S3). This negative lag correlation goes against the preferred maintenance of the AO sign between winter and summer (Ogi et al., 2004). Determining the sign of extreme erosion rates, solely based on the information of switches of the AO sign between winter and summer, yields accuracy rates of 70% for Bykovsky, 67% for Muostakh-N, and 33% for Muostakh-NE. Breaking down by sign of rate anomalies, true positive rates (sensitivity) are 80%, 75%, and 25%, while true negative rates (specificity) are 60%, 60%, and 38% for Bykovsky, Muostakh-N, and Muostakh-NE, respectively. Therefore, the linkage from the switch in sign between winter and summer AO to coastal erosion rates is more clearly noticeable at Bykovsky and Muostakh-N, than at Muostakh-NE.
Differences in results observed between Bykovsky and Muostakh may be due to local differences in geomorphological properties, which determine the prevailing erosion mechanism at each site. Thermo-denudation (TD) is identified as subaerial erosion primarily driven by temperature changes, hence ground-ice melt, permafrost thaw, and the consequent ground destabilization. Thermo-abrasion (TA) is characterized by the action of the kinetic energy of waves at undercutting coastal cliffs, leading to formation of notches at the land-sea interface, and subsequent rupture and fall of entire coastal blocks onto the shallow surf zone, which are then removed by ocean currents and waves. In TA-dominated coasts, a shoreline change is only apparent at the event of a block rupture, which may suddenly result after several months of mechanical abrasion of the cliff base by waves. Therefore, the TA component may add a delay to the shoreline position in response to external drivers. Although the relative role of TD and TA often shows high variability within small areas, Günther et al. (2015) showed that both processes were similarly important to total coastal erosion rates at Muostakh on average, and at Muostakh Cape specifically, in the period between 2010 and 2013. Lantuit et al. (2011) showed that different backshore landforms at Bykovsky are significantly associated with erosion rates. They described larger means and variability at coastal stretches affected by alases and retrogressive thaw slumps, which is the case of our key monitoring site at Bykovsky, at the southeastern portion of the Peninsula, facing the Laptev Sea to northeast. In this specific coastal segment, TD may play a more important role than TA in determining the total coastal erosion mean annual rates.
Both monitoring sites present ice-rich morphology and are thus prone to thermally driven erosion. The volumetric ground-ice content has been estimated at 79% for Bykovsky Peninsula by Fuchs et al. (2018) and at 87% for Muostakh Island by Günther et al. (2015). From the ACD database (Lantuit et al., 2012), the ground-ice content is estimated at 60% at both locations, larger than both the Laptev Sea and Arctic coast means of ∼23% and ∼13%, respectively ( Figure S5). Each of the two sites may also experience locally different weather conditions. Therefore, differences in response to the proposed drivers are also expected among sites.
Not only does coastal erosion respond to climate change, but also do sea ice and the AO. Pavlidis et al. (2007) estimated a 1.5-fold to 2.5-fold increase of Arctic coastal erosion rates by 2100, especially pronounced in areas where the largest changes in sea ice are expected. Barnhart et al. (2016) showed that the duration of the ice-free season at the Laptev Sea coast is projected to leave the normal range of variability by the end of the 21st century. Cai et al. (2018) showed a projected positive trend in the summer AO index between 2030 and 2100. Both studies followed a high-emission scenario (RCP 8.5, IPCC, 2013). Therefore, the relative role of the here depicted main modes of coastal erosion variability may change in the future with climate.

Conclusion
We identify the main modes of coastal erosion variability at the Laptev Sea based on almost 30 years of coastal erosion in-situ observations from Bykovsky Peninsula and Muostakh Island. The first mode accounts for 56.5±16.0% of the total variance and seems to respond to decadal-scale changes in winter SIC, predominant over much of the eastern marginal Arctic Seas, including the Laptev and East Siberian Seas. The second and higher-frequency mode accounts for 28.6 ± 8.1% of the total variance and can be associated with the AO during winter and summer. Together, the first two modes explain 85.1 ± 24.1% of the total variance; thus, most of the observed variability may rely on large-scale drivers.
SIC averaged over the Laptev Sea during winter (FMA) is an indicator for sea-ice melt start and end of freeze-up dates and hence modulates the duration of the OWS between two summer measurements. This is the first-order driver of coastal erosion variability at Bykovsky and Muostakh. The FMA AO likely contributes to coastal erosion by anticipating the onset of the melt season. Negative JJA AO conditions are associated with decreasing summer SIC anomalies, allowing increased fetch and wave activity, and with increasing surface warming, contributing to ground-ice melt. Therefore, the JJA AO would play a role on increasing Arctic coastal erosion by setting large-scale conditions favorable to both TA and TD. We also show that extreme erosion years follow a switch in sign from AO+ in FMA to AO− in JJA and vice versa. Taking only these years, significant negative correlations between the FMA AO and JJA AO indices emerge, indicating a combined effect of the two. We suggest that differences in the contribution of the proposed drivers between sites are attributed to local characteristics, such as ground-ice content, backshore material, and prevailing coastal erosion mechanism.
Our findings stand as an initial step towards bridging Arctic coastal erosion and large-scale climate variability. This is appealing because modern ESMs cannot physically represent coastal erosion at its small scale, but do represent synoptic-to planetary-scale mechanisms, such as pressure and wind changes, to which our proposed drivers directly respond. Also, the large-scale first-order drivers of coastal erosion may be common to other regions in the Arctic. Thereby, we recommend that Arctic coastal erosion start being considered in historical and future climate projections.