Increase in Arctic coastal erosion and its sensitivity to warming in the twenty-first century

Arctic coastal erosion damages infrastructure, threatens coastal communities and releases organic carbon from permafrost. However, the magnitude, timing and sensitivity of coastal erosion increase to global warming remain unknown. Here we project the Arctic-mean erosion rate to increase and very likely exceed its historical range of variability before the end of the century in a wide range of emission scenarios. The sensitivity of erosion to warming roughly doubles, reaching 0.4–0.8 m yr−1 °C−1 and 2.3–4.2 TgC yr−1 °C−1 by the end of the century. We develop a simplified semi-empirical model to produce twenty-first-century pan-Arctic coastal erosion rate projections. Our results will inform policymakers on coastal conservation and socioeconomic planning, and organic carbon flux projections lay out the path for future work to investigate the impact of Arctic coastal erosion on the changing Arctic Ocean, its role as a global carbon sink, and the permafrost–carbon feedback. Coastal erosion in the Arctic is caused by permafrost thaw and wave abrasion enhanced by sea ice melt, both of which will increase under climate change. Projections of erosion rate across the Arctic indicate that mean erosion rates will rise beyond historical precedent over the ﻿twenty-first century.

A rctic coastal erosion is caused by a combination of thermal and mechanical drivers. Permafrost thaw and ground-ice melt lead to soil decohesion and slumping, while surface ocean waves mechanically abrade the Arctic coast 1 . Sea-ice loss expands the fetch for waves 2,3 and prolongs the open-water season, increasing the vulnerability of the Arctic coast to erosion 4,5 . In the past decades, coastal retreat rates have increased throughout the Arctic, often by a factor of two or more [6][7][8][9][10] . The historical acceleration of erosion in the Arctic is linked with the observed decreasing sea-ice cover 2,4,11 , and increasing air surface 12,13 and permafrost temperatures 14 . As for the future, Arctic surface air temperature is projected to exceed its natural range of variability within the next decades 15 . Arctic sea ice decline has already exceeded natural variability 15 , and summer ice-free conditions are projected by mid-twenty-first century 16 . New regimes of surface waves are also projected in the Arctic Ocean and along the coast [17][18][19] . Consequently, Arctic coastal erosion rates are expected to increase in the coming decades. However, the extent of this increase is still unknown, as no projections of Arctic coastal erosion rates across the pan-Arctic scale are available. To fill this gap, we present twenty-first-century projections of coastal erosion at the pan-Arctic scale.
The thawing of permafrost globally releases organic carbon (OC) and increases atmospheric and oceanic greenhouse gas concentrations, feeding back to further warming [20][21][22][23] . Arctic coastal erosion alone releases about as much OC as all the Arctic rivers combined 23,24 , fuelling about one-fifth of Arctic marine primary production 25 . Despite consistent improvements in the representation of permafrost dynamics 26,27 , the current generation of Earth system models (ESMs) does not account for abrupt permafrost thaw, which may cause projections of OC losses to be largely underestimated 28,29 . Arctic coastal erosion is one form of abrupt permafrost thaw 22 and is a relevant component of the Arctic carbon cycle 23,30 . Nonetheless, it has not been considered in climate projections so far. The scale mismatch between Arctic coastal erosion and modern ESMs requires the development of holistic models that account for the key large-scale processes to bridge this gap [30][31][32] .
In this study, we present a novel approach to represent Arctic coastal erosion at the scales of modern ESMs. We develop a semi-empirical Arctic coastal erosion model combining observations from the Arctic Coastal Dynamics (ACD) database 33 , climate reanalyses, ESM and ocean surface wave simulations (see Methods for details). Our model considers the main thermal and mechanical drivers of erosion as dynamical variables, represented by yearly-accumulated positive temperatures and significant wave heights, and constant ground-ice content from observations. Our approach allows us to make twenty-first-century projections of coastal erosion at the pan-Arctic scale. We quantify the magnitude, timing and sensitivity of Arctic coastal erosion and its associated OC loss in the context of climate change.

Emergence of Arctic coastal erosion
We project the Arctic-mean coastal erosion rate to increase from 0.9 ± 0.4 m yr −1 during the historical period (1850-1950) to 1.6 ± 0.5, 2.0 ± 0.7 and 2.6 ± 0.8 m yr −1 by the end of the twenty-first century (2081-2100) in the context of anthropogenic climate change according to the socioeconomic pathway (SSP) scenarios SSP1-2.6, SSP2-4.5 and SSP5-8.5, respectively (Fig. 1a). This translates to an increase in the Arctic-mean coastal erosion rate by a factor of between 1.8 and 2.9 by the end of the century with respect to the historical period. The SSP2-4.5 and SSP5-8.5 scenarios describe moderate and high challenges to mitigation and adaptation to climate change 34 and, consequently, also moderate and high increases in radiative forcings due to greenhouse gas emissions 35 , respectively. The current cumulative CO 2 emissions is encompassed in the SSP2-4.5 and SSP5-8.5 range 36 . Scenario SSP1-2.6 follows a narrative of sustainability, where development is directed to environmentally friendly technologies, and society gradually transitions to renewable energy sources facilitated by international cooperation 34 . In the SSP2-4.5 and SSP5-8.5 scenarios, the Arctic-mean erosion rate continues increasing in the second half of the century, while in scenario SSP1-2.6, it stabilizes and starts showing decreasing trends.
We find it likely (≥66% probability) that the Arctic-mean erosion will exceed its historical range by around 2023 in all scenarios, and very likely (≥90% probability) by between 2049 (SSP5-8.5) and 2073 (SSP1-2.6), considering the largest distributions of uncertainties in our projections (that is, ensemble spread and erosion model uncertainties; Fig. 1b). The emergence of the Arctic-mean erosion rate from the historical range would very likely have happened by around 2010, if we take only the ensemble spread to define the historical range. Significant differences in projections between scenarios are only noticeable in the second half of the century, after a probable emergence from the historical range. Our erosion time-of-emergence estimates reflect those of its drivers, which take place around mid-twenty-first century (Fig. 1c,d), in accordance with previous studies 15,16 .
Arctic coastal erosion, often in combination with other climaterelated coastal hazards (that is, sea-level rise, permafrost thaw and storms), causes relevant socioeconomic damage 37,38 , including irreversible heritage data loss 39,40 . In Alaska, entire villages are already facing the need for relocation 41,42 . By the end of the century, different emission scenarios lead to substantially different projections in terms of societal impacts. The Arctic-mean erosion rate in SSP5-8.5 could reach values about twice as high as those in SSP1-2.6 by 2100, especially given the opposing projected trends. The Arctic-mean coastal erosion continues to steadily increase by the end of the century in SSP5-8.5, while it stabilizes in SSP1-2.6. Such differences in pace are decisive for the timing of adaptation from coastal communities.

Spatial variability of erosion
The thermal and mechanical drivers of erosion explain about 36-47% of its observed spatial variability in multiple linear regression models. On the one hand, wave exposure, combined with ground-ice content, best explains the spatial variability of erosion in most of the coastal segments (r = 0.69 ± 0.12, mean ± 2σ; Fig. 2b), where erosion is not extremely high ( ~ 90th percentile, < 2.5 m yr −1 ). The local wave exposure information indeed integrates important sources of erosion variability. Not only does wave exposure promote cliff abrasion and subsequent sediment transport, but it is also proportional to open-water season duration, which has been suggested to be the first-order driver of coastal erosion rate variability 2,32 . In addition, sea-ice melt, and thus increasing open-water season duration, responds to increasing surface air temperature, which also drives permafrost thaw and thus erosion by thermo-denudation. On the other hand, spatial differences among segments of extremely high long-term erosion rates are best characterized by thawing temperature exposure combined with ground-ice content (r = 0.61 ± 0.42; Fig. 2c). This suggests that thermo-denudation plays a more important role in driving coastal erosion rates at extreme-erosion segments than at non-extreme ones. Among both extreme and non-extreme erosion segments, ground ice adds explanatory power, as it increases the susceptibility of permafrost to thaw and hence erosion. Our results are in accordance with previous work, which reported spatial correlations between ground-ice content and erosion rates 33 . Strong temporal correlations between erosion and thawing temperature exposure have also been reported in the Mackenzie Delta region of the Canadian Beaufort Sea 43 , and in Muostakh Island of the Laptev Sea 8 , where erosion rates often range between 10 and 20 m yr −1 (refs. 8,44 ). We further combine the temporal evolution of the Arctic-mean erosion with its spatial distribution to make spatially resolved projections of erosion (Fig. 2d).
The geographical distribution of low and high-erosion segments does not change substantially from observations over time in our projections, which is partially a consequence of our model design,   as explained by the three following reasons. First, we assume that the spatial model coefficients, empirically determined, remain unchanged throughout our simulations. Second, ground-ice content, an explanatory variable in our regression model, is also assumed constant over time. Third, our regression model accounts for only a fraction of the spatial variability in erosion, and may thus underestimate larger spatial changes that may occur over time.
Moreover, and independent of model design, local anomalies of the dynamical variables (that is, local wave and thawing temperature exposure) are smaller in magnitude than their Arctic-mean increase. Therefore, our modelled changes in the spatial variability of erosion are small in comparison with its Arctic-mean increase. Nonetheless, our modelled spatial spread of erosion increases with time (Fig. 2e). The spread of erosion rate distributions increases from 3.6 m yr −1 (0-3.6, 5th-95th percentile range) in the historical period to 3.9 (0.9-4.8) and 4.2 (1.4-5.7) m yr −1 in the SSP2-4.5 and SSP5-8.5 scenarios, respectively. In scenario SSP1-2.6, the distribution spread does not increase, but shifts with respect to the historical period (3.6 (0.9-4.5) m yr −1 ). Temporally resolved erosion rate observations are rare, often sparse in time, and only available at a relatively small number of locations 10 . Only with such observations, temporally resolved and at the pan-Arctic scale, would empirical models be able to better constrain the temporal evolution of spatial variability of coastal erosion.

Spatial variability of organic carbon losses
The pan-Arctic OC loss from coastal erosion increases from 6.9  Table 2). For the present-day climate (that is, the period for which erosion observations are available), we estimate a pan-Arctic OC loss from coastal erosion of 8.5 (3.3-13.7) TgC yr −1 . Both our simulated present-climate mean and uncertainty range are comparable with previous estimates from observations 24,33 . Our projections suggest that pan-Arctic OC flux will increase by a factor of between 1.3 and 2.0 with respect to the present-day climate, or by a factor of between 1.7 and 2.5 by 2100 with respect to the historical period. SSP2-4.5 yields an increase by a factor of about 1.5-1.6 by the end of the century. The Laptev and East Siberian Seas (LESS, Fig. 2a) together account for about three-quarters of the pan-Arctic OC losses in our simulations, in accordance with observations-based estimates 24 . This also holds true for future scenarios. The reason for the relatively high OC fluxes from the LESS coast is twofold. First, the region comprises coastal segments of extremely rapid erosion, often between 10 and 20 m yr −1 (refs. 8,44 ). Second, the LESS coast is dominated by Yedoma ice-complex deposits, where ground-ice concentration reaches more than 80% of soil volume 8,45 , and organic-carbon content is extremely high, reaching about 5% of weight 33 . The role of coastal erosion in driving sediment fluxes from the LESS coast is relatively large in comparison with other Arctic marginal seas. Coastal erosion is estimated to contribute about twice as much as riverine discharge to the total sediment flux from the Laptev Sea coast 46 . Conversely, in the Canadian Beaufort Sea, the riverine contribution-modulated by the Mackenzie River-would be about 10 times larger than that from coastal erosion. More recent work 24 qualitatively confirmed the early estimates, and further suggests a modern increase in sediment and OC fluxes from coastal erosion in comparison with pre-industrial values across the Arctic. From the LESS, we simulate a present-climate OC flux of 6.5 (2.4-10.6) TgC yr −1 , comparable to the 2.9-11.0 TgC yr −1 range estimated by Wegner et al. 24 , and encompassing the ACD value of 7.7 TgC yr −1 . In an extensive investigation over the LESS continental shelf, Vonk et al. (2012) 23 determined that about 20 TgC yr −1 are buried in the LESS sediment, which would originate from a combination of coastal and seafloor erosion. Accounting for degradation before burial and assuming an equal contribution from coastal and subsea erosion, about 11 (7-15) TgC yr −1 would be released by coastal erosion alone. The LESS estimate of Vonk et al. (2012) 23 is 43-57% larger than other observations-based estimates 24 and about 69% larger than our present-climate modelled value. These differences are probably due to extensive and high-resolution sampling, allowing for more accurate upscaling 23 . However, the uncertainties associated with the contribution of coastal and subsea erosion encompass our modelled range (their Supplementary Table 6 23 ). Therefore, an underestimation from our side is not conclusive. Coastal erosion is estimated to sustain about one-fifth of the total Arctic marine primary production at present-climate conditions 25 . Therefore, the projected additional OC loss could have a substantial impact on the Arctic marine biogeochemistry. However, the fate of the organic carbon released by Arctic coastal erosion is currently under active debate. Field work has shown that between about 13% and 65% of the OC released into the ocean by coastal erosion could settle in the marine sediment [48][49][50] , slowing down remineralization. In the sediment, organic matter degradation would then take place at millennial time scales 51 . However, in the shallow nearshore zone, resuspension driven by waves and storm activity increases the residence time of OC in the water column, and allows for more effective remineralization 52 . Moreover, partial degradation of the eroded material takes place before it enters the ocean, releasing greenhouse gases directly to the atmosphere 22,23,53 . The OC degradation timescale thus also depends on its transit time onshore 53 . It is therefore challenging to determine short-term impacts from the projected additional OC fluxes from coastal erosion, as large uncertainties still remain regarding pathways of OC degradation.
We partition the uncertainty sources in our projections between three sources: ensemble spread, temporal and spatial erosion model components (see Methods). Our erosion model contributes the most to the uncertainties in our simulations: from about 76% of the total uncertainty range in the historical period and up to 97% by the end of the century in SSP5-8.5. The ensemble spread is responsible for the remaining 24% of the total uncertainty during the historical period, and for only 3-6% of the total range at the end of the future Modern climate (1950-2010) Lantuit et al. 33 Wegner et al. 24 Vonk et al. 23 Historical (1850-1900) SSP1-2.6 (2081-2100) SSP2-4.5 (2081-2100) SSP5-8.5 (2081-2100)  scenarios. The spatial component of the erosion model accounts for about half of the total range of uncertainties on average, without significant changes in proportion over time. The fraction of uncertainties stemming from the temporal model component increases from about 33% of the total range in the historical period to about 55% by the end of the century in SSP5-8.5 due to the increasing magnitude of the erosion drivers. The distribution of sources of uncertainties in our projections is qualitatively similar between the pan-Arctic and the regional totals.

Sensitivity of erosion and carbon losses to climate change
The sensitivity of Arctic coastal erosion to climate change increases with time in our simulations, and is tightly related with the Arctic amplification (AA) 12 after its onset. Arctic coastal erosion increases more rapidly in response to increasing global mean surface air temperature (SAT) in the future scenarios than it does in the historical period. Before the mid-1970s, neither global nor Arctic mean SAT decadal trends are consistently significantly positive (Fig. 4a). During this period, the correlation between the Arctic-mean erosion rate and the Arctic mean SAT is weak (r = 0.26 ± 0.29, mean ± 2σ range; Fig. 4b). However, after the 1970s, correlations between erosion and Arctic SAT increase substantially (SSP1-2.6: r = 0.73 ± 0.09; SSP2-4.5: r = 0.68 ± 0.18; SSP5-8.5: r = 0.93 ± 0.06, 2081-2100 means), driven by the concurrent increasing trends. This turning point is also marked by the AA onset, when the Arctic SAT starts increasing at a faster pace than the global SAT, that is, the AA factor is consistently significantly larger than 1 (Fig. 4c), except in the second half of the century in the SSP1-2.6 scenario, when global and Arctic SAT trends approach zero. Therefore, the sensitivity of erosion to global SAT reflects the sensitivity of Arctic SAT to global SAT-quantified as the AA factor-after the AA onset, given the strong correspondence between erosion and the Arctic SAT at that time. The sharp increase in erosion sensitivity and the AA factor to their maximum values in the early 2000s is a signature from the so-called 'hiatus' in global warming 54 . Global mean SAT stalls between the late 1990s and the early 2010s, while the erosion drivers continue to increase. Sensitivity values level off in the second half of the twenty-first century, when global mean SAT trends Running-window lengths are 20 yr in all plots, except in f, which has yearly resolution and no filtering. Different window lengths show qualitatively similar results (not shown). The AA onset (dashed pink line) takes place in 1976, when the Arctic SAT increases at a faster pace than the global mean SAT, that is, the AA factor is larger than 1. After the 1970s, the AA factor is consistently significantly larger than 1, except for late twenty-first century in the SSP2-4.5 and SSP1-2.6 scenarios, when global and Arctic mean SATs decelerated and 20 yr trends are momentarily similar.
decelerate. End-of-century sensitivities are lowest in the SSP2-4.5 scenario, when Arctic SAT trends decrease sharply to reach the also consistently decreasing global SAT trends, and the AA factor approaches one. To avoid the effect of the warming hiatus, we quantify erosion sensitivity considering the historical period until before the AA onset, and during the last 50 yr in the scenario simulations.
The sensitivity of the Arctic-mean erosion rate to global mean SAT increases significantly from 0.18 ± 0.31 m yr −1 °C −1 on average during the historical period until 1975, to at least double to between 0.40 ± 0.16 and 0.48 ± 0.21 m yr −1 °C −1 during the second half of the twenty-first century following the SSP2-4.5 and SSP5-8.5 scenarios, respectively. This translates to an increase in the sensitivity of OC losses to climate warming from 1.4 TgC yr −1 °C −1 in the historical period until 1975 on average, to between 2.3 and 2.8 TgC yr −1 °C −1 following the SSP2-4.5 and SSP5-8.5 scenarios, respectively. In scenario SSP1-2.6, the increase in sensitivity is steeper, reaching 0.79 ± 0.57 m yr −1 °C −1 , or about 4.2 TgC yr −1 °C −1 by the end of the century. However, all sensitivity estimates in SSP1-2.6 are associated with relatively large uncertainties, especially in the second half of the century, as global SAT trends approach zero. In addition, the AA factor decreases in magnitude to non-significant values and thus plays a smaller role in driving erosion in SSP1-2.6. The sensitivities of both the thermal and mechanical drivers of erosion increase at the AA onset (Fig. 4d,e), especially in scenarios SSP2-4.5 and SSP5-8.5, where they remain significantly positive until the end of the century, as opposed to the historical period. The mid-century peaks in the sensitivity of the erosion drivers, associated with large uncertainties, are a consequence of the SAT trends crossing the zero line.
The increase in the Arctic-mean erosion is larger in scenarios where erosion is mainly driven by thermal abrasion (TA), in comparison with thermal denudation (TD)-dominated scenarios (see Box 1). Sea-level rise increases the vulnerability of the Arctic coast to TA by increasing the inundated fraction of coastal cliffs, making them more susceptible to abrasion by ocean surface waves, niche formation and block failure. Similarly, but in a more episodic timescale, storms also increase coastal erosion by TA due to water level rise caused by surface winds and ocean surface waves. Attributing 90% of Arctic-mean erosion to mechanical contribution, we depict a scenario in which erosion is mainly driven by TA (q = 0.9 in Supplementary Table 1). This scenario yields Arctic-mean coastal erosion rates about 10% (SSP1-2.6) and (SSP5-8.5) 20% larger by 2100 in comparison with a scenario in which the role of TA is as large as that of TD and does not increase with time (q = 0.5). Increasing the weight attributed to TA necessarily increases erosion, as TA increases more than TD with respect to the historical period (Fig. 4f). Our wide range of proportionality factors should thus account for the potential increase in erosion stemming from an increasing role of TA with time by the end of the century. The increasing role of TA could thus be attributed to a combination of factors, including increasing sea level, or the increasing sensitivity of erosion to wave abrasion and storms. The effect of storms, both episodic effects and long-term changes, are directly considered in our model in the mechanical driver of erosion, which consists of yearly-aggregated daily wave information. Wave heights are forced by surface winds, which also increase the water level at the coast during storms which, in turn, increase the vulnerability of the coast to TA.
The sensitivity parameters are useful tools for assessing the state of Arctic coastal erosion increase and the associated OC fluxes at intermediate states or policy-based targets of global warming. It must be noted, however, that the sensitivity parameters assume linear relationships between the forcing and outcome variables 55 . Similarly, our model does not account for nonlinear, or synergistic effects that could emerge from the combination of the thermal and mechanical drivers of erosion. One could thus interpret our projections as a lower-bound estimate that accounts for the additive effect of its main drivers. We assume that a linear combination of the two main drivers of erosion provides us with a first-order large-scale estimate of the temporal evolution of Arctic coastal erosion, associated with a clear range of uncertainties.
Our erosion model is relatively simple compared with higher-resolution and process-based strategies (for example, refs. [56][57][58][59][60][61] ; see supplementary file for a summary of previous modelling work), and does not intend to represent all processes, often of fine spatial scale (order of metres or less), associated with the erosion of the Arctic coast. High-resolution process-based modelling strategies are useful for predictions at the local scale (for example, ref. 62 ), or for investigating individual aspects of erosion in isolation (for example, ref. 59 ). However, it is unfeasible to run such detailed models for climate projections at the pan-Arctic scale. Here we empirically parameterize coastal erosion at scales compatible with the resolution and mechanisms represented in ESMs (order of tens or hundreds of kilometres). Despite the limitations at reproducing erosion at the local scale, our semi-empirical approach allows us to make pan-Arctic projections of coastal erosion and its associated OC fluxes, and thus make first-order estimates of the magnitude, timing and sensitivity of their increase to global warming. We thereby highlight the need for further work to constrain our relatively large uncertainties.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41558-022-01281-0.

Box 1 | Thermal and mechanical drivers
Arctic coastal erosion is typically caused by a combination of TD and TA 1 , which act together to thaw permafrost, melt ground ice, abrade and transport coastal material offshore. We take yearly-accumulated daily positive temperatures and significant wave heights to represent TD and TA: hereafter, the 'thermal' and 'mechanical' drivers of erosion, respectively. As various landform types compose the Arctic coast, the relative contributions of the thermal and mechanical drivers differ at the local scale. Erosion is predominantly thermally driven at retrogressive thaw slumps, observed at the Bykovsky Peninsula, Laptev Sea 64 and in the Mackenzie Delta region (Beaufort Sea) 43,65,66 , for example (Fig. 2a), as wave abrasion and sediment transport from ocean currents play a secondary role in coastal retreat in such formations. Erosion is also predominantly thermally driven in enclosed bays and in coastal segments protected by spits and barrier islands, where the fetch of ocean waves is limited 67 , although barrier islands themselves are often susceptible to wave abrasion 68 . In contrast, erosion of ice-rich coastal cliffs occurring extensively along the Beaufort and Laptev Sea coasts, for example [6][7][8] , requires the mechanical action from ocean waves to open notches at the land-sea interface, causing the subsequent collapse of often still-frozen large blocks of permafrost. In some locations, the relative contributions of the thermal and mechanical drivers are more balanced than described above. At Muostakh Island in the Laptev Sea, for example, TD and TA are estimated to similarly contribute to maintain erosion rates of up to 25 m yr −1 (refs. 8,44 ).

Data.
Arctic coastal observations. We used the ACD database 33 as our observational reference. The ACD compiles several sources of data and provides a list of variables for a total of 1,314 coastal segments along the Arctic coast, including: long-term mean erosion rates, organic carbon concentration, soil bulk density, ground-ice fraction, mean elevation and length. From the 1,314 segments, we took the 306 classified as erosive and non-lithified, which excludes segments from the rocky coasts of Greenland and the Canadian archipelago, and other segments that present stable or aggrading dynamics. We also selected segments containing excess ice, which excludes all the non-erosive segments from Svalbard, for example. We used a subset of 72 coastal segments, classified with a 'high-quality' flag in the ACD, to train our semi-empirical model. The model was validated and then forced over the entire set of 306 coastal segments for the pan-Arctic simulations.
Reanalysis. We took 2 m air temperature, significant wave height and sea-ice concentration data from ERA20C reanalysis 69 as empirical variables to train and validate our coastal erosion model. Data were taken for the same periods for which the erosion rates are provided in the ACD. Data have ~1.12° (temperature and sea ice) and 1.5° (waves) horizontal resolution. We assigned the closest land grid cell in ERA20C from its atmospheric grid to ACD segments, and two rows of adjacent cells from the ocean grid.
Climate projections. We used a 10-member ensemble of simulations from the Max Planck Institute Earth System Model (MPI-ESM) version 1.2 in its low-resolution configuration 70 performed for the Coupled Model Intercomparison Project phase 6 (CMIP6) 35 . In this configuration, the atmospheric component ECHAM6 71 works on the T63 setup, which has 1.875° of nominal horizontal resolution on a Gaussian grid, and provides about 100 km grid spacing along the Arctic coast. The ocean component MPIOM 72 works on the GR15 setup, which consists of a curvilinear grid of 1.5° of nominal resolution. This corresponds to a horizontal grid spacing ranging from 15 km around Greenland to about 150 km in the tropics, and about 80 km along the Arctic coast. We used the historical simulations (1850-2014) and three future shared socioeconomic pathways (SSPs) for the twenty-first century projections (2015-2100), namely: the SSP1-2.6, SSP2-4.5 and the SSP5-8.5, which represent a wide range between low-and high-end emission scenarios 34 . This range of scenarios is realistic in terms of current cumulative CO 2 emissions 36 . From the MPI-ESM, we took 2 m temperatures and sea-ice concentrations as input for our erosion model. Surface winds and sea-ice concentration data were also used to force ocean surface wave simulations.
Ocean wave simulations. We used the wave model WAM 73 to generate a 10-member ensemble of global waves for historical, SSP1-2.6, SSP2-4.5 and SSP5-8.5 scenarios, forced by the MPI-ESM ensemble. In our setup, WAM has 1° grid resolution and is forced with daily sea-ice concentration (threshold of 15% to define open-water), 6-hourly 10 m winds, and a realistic ETOPO2-based bathymetry as boundary conditions. From WAM simulations, we took daily-mean significant wave heights (H s (m)), defined as the mean over the highest one-third of the total wave-height distribution.
Semi-empirical Arctic coastal erosion model. We present a simplified model for Arctic coastal erosion, compatible with the scales of Earth system models. Our model considers the dominant physical thermal and mechanical drivers of erosion, also referred to as TA and TD 1 . The model is constrained to only simulate erosion in the presence of ground ice and in the absence of coastal sea ice. We used an empirical approach to quantify the relationship between the physical drivers, constraints and erosion rates, by comparing the observations from the ACD with ERA20C reanalysis. The empirically estimated parameters were then applied to all coastal segments, which provides us with erosion rates in the pan-Arctic scale. Our model has yearly time resolution, and the spatial resolution follows the definitions of the ACD coastal segments.
The total erosion E(x,t) (m yr −1 ), defined in every year t and coastal segment x(lat, lon), is given as a combination of a temporal and a spatial component.
The temporal component represents the temporal evolution of the Arctic-mean erosion E(t) (m yr −1 ). The spatial component ΔE(x,t) (m yr −1 ) represents local departures from the Arctic mean at every year and coastal segment, providing spatially distributed values of erosion. Hereafter, we use ' Arctic mean' , denoted by the overline, to refer to means along the Arctic coast. All data associated with ACD coastal segments were weighted by segment lengths in the computation of means.
The temporal component. The temporal component of our model is a linear combination of Arctic means of the thermal and mechanical drivers of erosion.
The thermal driver of erosion is represented by Arctic-mean yearly-accumulated daily-mean positive 2 m air temperatures T(t) (°C d yr −1 ), also commonly known as positive degree-days or thawing degree-days. The mechanical driver of erosion is represented by Arctic-mean yearly-accumulated daily significant wave heights H(t) (m d yr −1 ).
We empirically estimated the linear coefficients a TA (m m −1 d −1 yr) and a TD (°C m −1 d −1 yr) by scaling the Arctic-mean physical drivers from ERA20C reanalysis, with the observed coastal erosion rates from the ACD. This was done for the reference time t obs , during which observations are available.
We assumed that the thermal and mechanical drivers a TD T(t) and a TA H(t) contribute in equal proportions to the Arctic-mean erosion during the reference time. We did this by setting the proportionality factor q to 0.5. We tested the sensitivity of our results to this assumption by making scenarios with q = 0.1 and q = 0.9. This sensitivity test yielded a maximum increase (decrease) of about 10-20% in projections to 2100 with respect to the equal-proportion case (Supplementary Table 1).
The spatial component. The spatial component of our erosion model calculated local erosion anomalies with respect to the Arctic-mean temporal evolution, and consisted of two multiple linear regression (MLR) models. We split the coastal segments in two groups by classifying them between 'extreme' and 'non-extreme' with respect to erosion, using 2.5 m yr −1 as a threshold (~90th percentile). We did not find a distinct separation between extreme and non-extreme segments in terms of geographical location (Fig. 2a), or coastal morphology. Both groups showed similar distributions of ground-ice content, mean cliff height, bathymetric profile, bulk density, as well as mean thermal and mechanical forcings derived from thawing temperature and ocean waves, for example (not shown). We tested a comprehensive number of combinations of dynamical and geomorphological parameters as explanatory variables in MLR models, simultaneously maximizing goodness-of-fit and penalizing model complexity (Supplementary Table 3). We fit MLR models using the usual ordinary least square method. The goodness-of-fit of models was assessed with the proportion of explained variance and root mean squared error (RMSE). Since increasing the number of combined explanatory variables necessarily increases the model fit and may lead to overfitting, we penalized model complexity by assessing the changes in the Akaike Information Criterion 74 (ΔAIC) in parallel. The best performing combination of covariates is the one which maximizes correlation (or proportion of explained variance) and minimizes RMSE and ΔAIC (Supplementary Fig. 2). We trained the spatial component of our erosion model only on those segments classified as 'high quality' with respect to erosion data. We included medium-quality segments to train the model for the high-erosion case to increase our sample size and thus also statistical robustness. We validated each combination of regression coefficients with unseen data by performing a leave-one-out cross-validation test. We used a Bootstrap approach with 10,000 sampling iterations to obtain distributions of model coefficients and hence their associated uncertainties.
Three variables composed the best performing combinations: (1) daily-mean thawing temperature exposure, expressed as the yearly-accumulated daily positive temperature divided by the number of positive-temperature days per year T day (°C yr −1 ), (2) daily-mean wave exposure, expressed as the yearly-accumulated daily significant wave heights divided by the number of open-water days per year H day (m yr −1 ), and (3) ground-ice content θ (% of soil volume). On one hand, combining ground-ice content with daily-mean wave exposure (θ+H day ) explained about 47% of the observed spatial variance among non-extreme (2.5 m yr −1 threshold) erosion segments (r = 0.69, 9-95th-percentile range: r = 0.60−0.78; Fig. 2b and Supplementary Fig. 3a). On the other hand, combining ground-ice content with the daily-mean thawing temperature exposure (θ+T day ) explained about 36% of the variance among extreme-erosion segments (r = 0.61, 9-95th-percentile range: r = 0.31−0.94; Fig. 2c and Supplementary Fig. 3a). The linear regression coefficients b obtained with the selected variable combinations were statistically significant (P < 0.01).
Swapping the combinations and groups, that is, using θ+H day for the extreme and θ+T day for the non-extreme erosion segments yielded overall poorer fits ( Supplementary Fig. 3a,b) and less robust estimation of regression coefficients ( Supplementary Fig. 3c-e). We also tested the sensitivity of these results to the choice of the threshold to define extreme erosion. Allowing for an overlap between the extreme and non-extreme segments by lowering the threshold to 2.0 m yr −1 , for example, increased the robustness of the T day regression coefficient estimate for the extreme group (Supplementary Fig. 3d) by increasing the number of data points, and yielded a similar fit to that of the higher threshold (θ+T day in Supplementary  Fig. 3a,b) and also similar ground-ice coefficients (θ+T day in Supplementary Fig. 3c).
Finally, the total erosion was constrained to the open-water period, and set to zero whenever and wherever sea-ice concentration (SIC) was above 15% at the coast. Combining the temporal (equation 2) and spatial (equation 5) components into our total erosion model (equation 1), conditioned by open-water and the extreme-erosion threshold, our model assumed the complete form: Bias correction. Before forcing the erosion model with MPI-ESM data, we adjusted the historical and scenario simulations for climate biases. The bias was removed between ERA20C data (used to estimate our model parameters) and MPI-ESM ensemble means at the coastal segments and reference periods from observations. The modelled distributions were shifted and scaled so that their means and spread fit those of ERA20C at the reference time.
Organic carbon fluxes. We translated linear erosion rates into volumetric erosion rates E vol (m 3 yr −1 ), sediment fluxes S (kg yr −1 ) and carbon fluxes C flux (kg yr −1 ), considering the mean geometry and ground properties of each coastal segment.
where L and h are the segments' mean length and elevation (m), θ is the ground-ice content (% volume), ρ is the soil bulk density (kg m −3 ) and C conc. is the organic carbon concentration (% weight). We integrated over the coastal segments: to obtain the total Arctic flux.
Sensitivity to climate change. We estimated the sensitivity of the organic carbon release by Arctic coastal erosion to climate change following the approach of Friedlingstein et al. 55 , but with a simplified set of tools. In their work, Friedlingstein et al. compared pairs of 'coupled' and 'uncoupled' simulations, where the increasing atmospheric CO 2 concentration either affected climate, or was neutral in terms of radiative effect. This pairwise comparison was necessary because the landatmosphere and ocean-atmosphere carbon fluxes respond to changes in both climate and atmospheric CO 2 concentrations. Therefore, the difference between their coupled and uncoupled simulations isolated the effect of the CO 2 -induced changes in climate on carbon fluxes from the effect of the changing atmospheric CO 2 concentration. In our case, changes in atmospheric CO 2 alone do not induce any Arctic coastal erosion response, if not by its radiative effect. An uncoupled simulation, where CO 2 does not induce a change in climate, would not yield any change in the organic carbon released by Arctic coastal erosion. Therefore, we can estimate the sensitivity of organic carbon release by Arctic coastal erosion to climate γ (TgC yr −1 °C −1) by comparing changes in global mean surface temperature and the resulting changes in carbon fluxes from erosion.
Probability and onset of emergence from the historical range. We defined the yearly probability density distribution of a modelled variable ψ as the normal distribution N(t) at year t. The mean of N(t) is the ensemble mean and its standard deviation is the ensemble standard deviation (plus the standard deviation of the distribution of erosion model uncertainties in specific situations, explained in the text). Similarly, the historical range of a modelled variable ψ is the normal distribution fitted to its average over the period 1850-1950, N hist . We calculated the area of distributions A hist = ∫N hist dψ and A(t) = ∫N(t)dψ to determine their overlap A hist ∩ A(t). We defined the probability of emergence from the historical range P(t), that is, the probability that N(t) will be different from N hist , as the fraction of A(t) that emerges from A hist : We defined the onset of emergence as the year when the ensemble mean is larger than μ + 2σ from historical range N hist .
Estimation of uncertainties. All ranges of uncertainties, except when clearly stated otherwise, were calculated with a Bootstrap method, which suits cases where the number of data is relatively small. From any vector X of arbitrary length, a large number (that is, 10,000) of vectors X i (i = 1, 2, … 10,000) was generated by sampling with replacement from X. The uncertainty of any statistics of X was estimated from the distribution of i realizations of the statistics obtained from X i .