Accelerated sea ice loss in the Wandel Sea points to a change in the Arctic’s Last Ice Area

The Arctic Ocean’s Wandel Sea is the easternmost sector of the Last Ice Area, where thick, old sea ice is expected to endure longer than elsewhere. Nevertheless, in August 2020 the area experienced record-low sea ice concentration. Here we use satellite data and sea ice model experiments to determine what caused this record sea ice minimum. In our simulations there was a multi-year sea-ice thinning trend due to climate change. Natural climate variability expressed as wind-forced ice advection and subsequent melt added to this trend. In spring 2020, the Wandel Sea had a mixture of both thin and—unusual for recent years—thick ice, but this thick ice was not sufficiently widespread to prevent the summer sea ice concentration minimum. With continued thinning, more frequent low summer sea ice events are expected. We suggest that the Last Ice Area, an important refuge for ice-dependent species, is less resilient to warming than previously thought. Arctic sea ice dynamics during August 2020 reveal that the Last Ice Area is more vulnerable to climate change and natural climate variability than previously thought, according to analyses of satellite data and sea ice model experiments.

I n August 2020, as the German Icebreaker Polarstern transited northward as part of the MOSAiC experiment, it followed a route across the Wandel Sea (WS) north of Greenland (Fig. 1a). This area is normally marked by compact, thick multiyear ice created by cold temperatures and onshore winds and currents 1 that stretches from the WS southwestward along the northern coast of the Canadian Arctic Archipelago. This area has been predicted to retain multiyear ice much longer than elsewhere 2,3 as the climate warms, and is therefore referred to as the "Last Ice Area" 4 (LIA, Fig. S1a). In fact, predictions of when sea ice-free summers in the Arctic will occur usually assume an area threshold of 1 × 10 6 km 2 , a number that implicitly accounts for the LIA 3,5 . The LIA is considered to be a last refuge for iceassociated Arctic marine mammals, such as polar bears (Ursus maritimus), ice-dependent seals such as ringed seals (Pusa hispida) and bearded seals (Erignathus barbatus), and walrus (Odobendus rosmarus) throughout the 21st century 6,7 . The LIA is also important for ivory gulls that breed in north Greenland 8 .
The first sign of change in the LIA emerged during the spring of 2018, when a large polynya formed in response to anomalous northward winds which drove sea ice away from the coast 9,10 . These winds were so strong that model simulations indicated that the polynya would have developed even with the thicker ice that had been present there several decades ago 11 . This serves as a reminder that not all notable events encountered in the Arctic are the result of climate change and much is yet to be discovered and understood.
But what about the summer of 2020? How unusual were the sea ice conditions that helped the Polarstern on its northward journey? Was this simply an example of interannual variability or did climate change contribute significantly to these ice conditions? And most importantly, is the summer of 2020 an early warning that the WS may not be as resilient to climate change as we expect? Here we use satellite data and sea ice model experiments to diagnose the mechanisms responsible for the summer 2020 record minimum and to put them into historical context. We also speculate about potential future changes in this area. We examine the interplay of weather-and climate change-driven changes in sea ice conditions and hypothesize about the potential short-and long-term implications for sea ice-associated species. Our focus is on the WS in the eastern LIA, although the implications of this study may be relevant to the entire LIA region.

Results and discussion
Observed summer sea ice conditions. The Polarstern's route was guided by satellite images showing extensive areas of open water and sea ice concentration (SIC) as low as 70% at 87N (Figs. 1a, S1b). We define our WS study area by 81.5°N-85°N, 10°W-50°W, the same area where we saw signs of change in February 2018 10 . Daily 2020 WS SIC drops below the 5th percentile of the 1979-2020 time series on July 25 and stays there almost until the end of August (Fig. 1b). August 14, 2020 constitutes a record low 52% SIC minimum (Fig. 1c). Several earlier years (e.g., 1985: 57%, 1990: 67%, and 1991: 62%) also show significant low SIC minima, although none as low as 2020.
Summer 2020 sea ice mass budget. Summer sea ice loss in any particular area happens in response to ice advection (i.e., dynamics) and ice melt (i.e., thermodynamics). Thus, to fully understand the causes of the summer 2020 WS sea ice loss event, we must examine these processes separately and in detail. To do this, we here use the Pan-arctic Ice Ocean Modeling and Assimilation System (PIOMAS; see "Methods"). For a time period Δt, a change in sea ice thickness (SIT) Δh is described by Δh/ Δt = F adv + F prod , where "ice advection" F adv is a short-hand term used here for SIT flux convergence. Positive (negative) ice advection indicates ice gain (loss). The quantity F prod is net ice production defined by F prod = F atm-ice + F bot where F atm-ice is surface ice melt due to atmospheric heating, and F bot is ice melt on the bottom and lateral floe surfaces due to ocean-ice heat flux (e.g. ref. 12 ).
First, we discuss ice loss from dynamics. Ice advection anomalies dominate monthly variability from 1979 to 2020, but with no long-term trend (Fig. 2a). Anomalous ice advection for June, July, and August (JJA) 2020 is substantially negative (i.e., SIT loss), although several earlier years also have anomalies of similar magnitude. July and August ice motion anomalies show ice being driven northwestward away from the WS (Fig. 2b). Daily F adv variability over the summer in 2020 (Fig. 2c) shows several strong negative events in early June, and mostly negative values for much of July through mid-August (see Fig. S2b for anomalies). Many, but not all of these F adv extrema correspond to peaks in surface stress as a result of strong wind events (Fig. S2a). The temporal evolution of cumulative SIT anomaly closely tracks the cumulative sea ice advection anomaly (Fig. 2d), underlining the importance of dynamics on SIT in the WS.
Next, we discuss ice loss from thermodynamics. In contrast to ice advection, Fig. 2a shows that 2020 JJA ice production is a record low value of −0.3 m/month (i.e., SIT loss). In fact, all years since 2016 show negative JJA ice production anomalies. Figure 3a shows that in 2020, F atm-ice is generally the largest contributor to ice melt until mid-July, when F bot becomes dominant as the ice recedes and more solar radiation enters the ocean. In order to better understand the physics of F bot , we write the ocean heat budget 12 for the upper 60 m (i.e., the layer that includes both the summer and winter mixed layer depths) of the WS region as: where ΔH is the change in ocean heat content over a time period Δt (here one day), F atm-ocn the net atmospheric heat flux at the ocean surface either in open water or under the ice cover (the latter referring to penetrating solar shortwave radiative flux), F bot the ocean-forced melting, and F ocdyn the sum of all oceanic dynamical heat fluxes at the lateral and bottom boundaries including convection and lateral and vertical advection and diffusion (see "Methods"). The dominant term in summertime F atm-ocn is solar shortwave radiation 13 which contributes to both bottom melt F bot and upper-ocean warming ΔH. Figure 3a shows how the terms in Eq. (1) vary over the summer of 2020 in the WS and how they interact to create the 2020 SIC anomaly. An increase in F atm-ocn starting in mid June increases ocean heat content H and results in warmer subsurface ocean temperatures that form a Near-Surface Temperature Maximum (NSTM) 12 (Fig. 3b). This heat mixes upward during periods of high surface stress, cooling the NSTM and melting ice (most notably during August 9-16 when the upper part of the NSTM above the temperature maximum loses heat). We note that F ocndyn is generally near zero, meaning that net ocean dynamics plays only a secondary role in ocean warming and bottom ice melt. This near-zero value is actually a sum of substantial warming via lateral heat flux convergence that is balanced by cooling via downwelling (see "Methods"). Total ice thickness during the time from July 25 and August 16th when we track the NSTM drops from about 3 m to 1.5 m (Fig. S3b) Fig. 2 Sea ice mass budget terms, ice thickness, and ice velocities. a Anomalies (relative to mean 1979-2020 June-August (JJA)) sea ice advection (F adv ) and ice production (F prod ) from 1979 to 2020 from the PIOMAS HIST run over the years from 1979 to 2020, b mean sea ice motion for JJA 2020 from the PIOMAS HIST run, c as in (a) but for daily JJA 2020 anomalies and including the sea ice thickness (SIT) anomaly, d as in (c) but in cumulative form. Budget terms F adv and F prod are given in units of meters of ice melt per month (panel a), meters of ice melt per day (panel c), and meters of ice melt (d).
It is important to recognize that ice advection and ice melt are not independent of each other. Negative advection (i.e., thickness divergence) creates thin ice and open water which then allows for more absorbed solar radiation, which in turn enhances melt: the classic ice-albedo feedback. This is evident in Fig. 3a and S2b, where F atm-ocn increases during periods of negative advection (e.g. during . This is a key result of our study: we find significant ice-albedo feedback in the WS, an area of usually thick, compact ice that traditionally should not sustain such a process 14 . The above analysis shows that the development of the record SIC anomaly in the WS was in large part driven by ice advection whose effect was amplified by anomalous melt. The fact that dynamics, specifically sea ice advection, in this area displays large interannual variability with no statistically significant trend suggests that low ice concentrations and low thickness during August 2020 may have been predominantly driven by dynamics. For such events, climate change attribution is often aided by separating the thermodynamic component which often more clearly shows a climate change signal 15 . In our case, we consider the long-term thinning of the ice cover to be dominated by thermodynamics 16 analogous to the sea surface temperatures identified in (ref. 15 ). Therefore, we next take a closer look at the sea ice conditions in spring 2020 (what we will call "initial conditions" for the events later that summer) relative to past years.
Initial conditions and the role of the ice thickness distribution. During spring 2020, ice accumulated in the WS (Fig. 4a, b) in response to anomalous advection (mostly in February; Fig. 4c, d).
As a result, ice thickness was near its 1979-2020 mean value by June 1 according to PIOMAS; Fig. 2c), and actually thicker than in recent years (2011-2019) as confirmed by the combined CryoSat-2/SMOS satellite product ( Fig. 4b) 12 . We might then wonder how this initial spring state evolved into the low SICs observed in August. As shown next, mean SIT does not tell the entire story: we must look at the sub-grid scale distribution of SIT which describes how much ice in each model grid cell is relatively thin or thick (see "Methods").
First, we consider how the SIT distribution changed over 2020, specifically from January through August 2020 (Fig. 5a). We see how thin ice in January transformed (via advection, Fig. 4c) to a distribution with significant thick ice while still retaining large amounts of thin ice. This substantial thick ice component is unusual in the context of recent years and is corroborated by observations that show large amounts of old ice in 2020 that are not evident in recent years (Fig. 5b).
Next, we consider the bigger picture i.e. longer-term changes in the SIT distribution over the period 1979-2020. During this time, the fraction of ice mass in PIOMAS's thicker ice categories decreased, while the fraction in the thin ice categories increased (Fig. 5c), in keeping with the general thinning trend observed in the Arctic 17,18 as well as in the LIA 10 . Thus, in recent years, thinner ice is the norm (Fig. 5d), which could lead to low summer SICs (such as seen in 2020) via reduced ice strength, increased susceptibility to advection, and enhanced shortwave absorption by the ocean and thus anomalous melt. In fact, one might ask, why haven't other recent years also shown low summer SICs?
The following section attempts to answer this question, looking The surface heat budget and upper ocean hydrographic profiles. a Upper (60 m) ocean heat budget terms for the WS over JJA 2020 from the PIOMAS HIST run, also with F adv at the top. ΔH/Δt is the change in upper ocean heat content over a time period Δt, F adv is sea ice advection, F atm-ocn is the net atmospheric heat flux at the ocean surface, F atm-ice is surface ice melt due to atmospheric heating, F bot is ice melt on bottom and lateral floe surfaces due to ocean-ice heat flux, F ocndyn is the sum of all heat fluxes due to ocean dynamics. Positive values denote warming. Surface stress (Nm −2 ) is calculated in the model at the top of the ocean, and so is a combination of air-ocean and ice-ocean stresses in partially ice-covered conditions. All quantities given in meters of ice melt per day (see "Methods"), b-e Temperature and salinity profiles from the PIOMAS HIST run for select dates during summer 2020. Shaded areas indicate the Near Surface Temperature Maximum, with maximum temperature (C) T max and ocean heat content (MJm −2 ) above the T max noted for the latter three dates.
specifically at the roles of springtime initial conditions vs. summer atmospheric forcing.
Model experiments: initial SIT vs. summer weather. Two model experiment ensembles (see "Methods" for further details) are analyzed in this section and compared to the historical simulation "HIST" previously described. First, we ask what role did the particular sea ice conditions found at the end of spring 2020 play in the extreme sea ice loss over the following summer? Experiment "INIT" addresses this question by starting 42 separate simulations using June 1 sea ice conditions from each of the years 1979-2020, while always using the atmospheric forcing only from 2020. We next ask, what role did the particular atmospheric forcings (wind, heat fluxes, etc.) found during summer 2020 play in the sea ice loss of that year? Experiment "ATMOS" addresses this question by starting 42 separate simulations all with June 1, 2020 sea ice conditions but then forcing the ice-ocean model with the atmospheric forcing from each summer of 1979-2020.
SIC for the HIST simulation is below the median of the INIT ensemble throughout the summer, dropping below the 10th percentile by the middle of June and staying there until the end of August (Fig. 6a). This shows that the spring initial conditions were key in producing the anomalous summer conditions in 2020. However, recent years with similar amounts of thin ice to 2020 (i.e., 2018 and 2019) also show low August SIC under 2020 atmospheric forcing. This means that if 2020 summer atmospheric conditions had occurred in 2018 or 2019, the WS would have seen similarly low ice concentrations as in 2020. It also suggests that the thicker and older ice present at the start of the summer of 2020 did little to prevent the record SIC anomaly in August.
The ATMOS ensemble (Fig. 6b) indicates that the impact of 2020 atmospheric forcing on that summer's SIC anomalies was minimal until late July, but then this forcing did indeed become important. HIST SIC drops below the 10th percentile of the ATMOS ensemble by the end of July, staying there through the end of August.
These ensemble experiments underline the importance of both spring sea ice and summer atmospheric forcing to August SIC. In summary, we find that: Spring ice conditions were mostly responsible for the summer SIC anomaly through the end of July, while the atmosphere was mainly responsible for driving SIC to a record low during August. Partitioning the impact of 2020 spring initial sea ice conditions vs. summer atmospheric forcing on the sea ice anomaly at the time of the WS sea ice minimum on August 14 (see "Methods") attributes~20% to the initial conditions while~80% is the due to the atmospheric forcing. Assuming that the increasing presence over the past 40 years of thin ice and open water at the beginning of the melt season (Fig. 5c) is primarily driven by climate change and that the summer atmospheric conditions were due to internal variability, the above 20/80 partitioning provides an approximate measure of the contributions of climate change and internal variability to the 2020 event (see further discussion in "Methods").
Is this 20% climate change signal significant or not? Extreme local events such as storms, heat waves, or floods are almost always dominated by dynamics driven by internal climate variability 15 . For example, flooding in New York City in response to Superstorm Sandy was on the order of 20% more extreme owing to long-term sea level rise 17,18 , a signal of the same magnitude we have detected in the present study. Using this example as a scale, we conclude that climate change was indeed a significant contributor to sea ice loss during the summer of 2020 in the WS.
Atmospheric variability in context. To put 2020 WS sea ice advection into a larger scale context, we consider here the fundamental modes of Arctic atmospheric variability, i.e., the Arctic Oscillation (AO), the Arctic Dipole Mode (ADA), and the Barents Oscillation (BO) 19,20 . Each of these correspond to the principal components (PCs) of empirical orthogonal functions computed from monthly mean sea level pressure fields north of 30°N (Figs. 7 and S4). During January-March 2020, when sea ice was advected into the Wandel Sea, sea level pressure over the Arctic was low, with a sea level pressure pattern similar to that found in 2017 when the Beaufort Gyre reversed 21,22 . The resulting onshore ice motion contributed to anomalously thick ice north of Greenland. At this time, the AO and ADA were both very high (the AO was in fact a record), a situation not found in any other year over the 41-year time series. Interestingly, summer 2020 conditions show the opposite, with ice motion westward away from the WS and the AO and ADA near record negative values. It seems clear that the anomalous 2020 WS wind forcing was associated with anomalous large-scale surface wind patterns.
Discussion and conclusions. While primarily driven by unusual weather, climate change in the form of thinning sea ice contributed significantly to the record low August 2020 SIC in the WS. Several advection events, some relatively early in the melt season, transported sea ice out of the region and allowed the accumulation of heat from the absorption of solar radiation in the ocean. This heat was mixed upward and contributed to rapid melt during high wind events, notably between August 9 and 16. Ocean-forced melting in this area that is traditionally covered by thick, compact ice is a key finding of this study. However, in some ways this should not be surprising given that this area (like most areas of the Arctic) has experienced a long-term thinning trend (Fig. 5c). Given the long-term thinning trend and strong interannual variability in atmospheric forcing, it seems reasonable to expect that summer sea ice conditions in the WS will likely become more variable in the future. In fact, mean SIT at the start of summer in 2018 and 2019 was even thinner than in 2020, which implies that with 2020-type atmospheric forcing, we might have seen even lower August SICs in those years, relative to that observed in 2020. In other words, SIC in the WS is now poised to plunge to low summer values, given the right atmospheric forcing. So, is the LIA in trouble?. The WS is a key part of the LIA, one that has recently experienced anomalous conditions 9,10 . We have shown that climate change-associated thinning ice in this region is a prerequisite for the record low ice concentrations seen in August 2020. Further, the unusually high SIT at the start of 2020 suggests that a temporary replenishment of sea ice from other parts of the Arctic may do little to protect this area from eventual sea ice loss. Recent work indicates that while western and eastern sectors of the LIA have distinct physics 11 , they both are experiencing long-term sea ice thinning and thus are both vulnerable to the processes discussed in this study. Our work suggests a re-examination of climate model simulations in this area, since most do not predict summer 2020-level low SICs until several decades or more into the future. Given that most climate models presently feature a subgrid-scale thickness distribution 23 its evolution over time in those models should be a focal point of future investigations (i.e., rather than simply focusing on grid-cell mean thickness). Coupled model simulations where atmosphere and ocean conditions are nudged to 2020 conditions would provide useful insights into the capabilities of the current generation of climate models to replicate our results. In addition, our results should be replicated with ice-ocean models using different resolutions and physics.
While the WS is only one part of the LIA, our results should give us pause when making assumptions about the persistence and resilience of summer sea ice in the LIA. Currently, little is known about marine mammal densities and biological productivity in the WS and the broader LIA. Recent studies indicate there may be some transient benefits for polar bears in areas transitioning from thick multi-year ice to thinner first year ice 24,25 , as biological productivity in the system increases 26 . However, this is largely the case in shallow water <300 m in depth and it is unclear if this will occur in multi-year ice regions elsewhere. The assumption that the LIA will be available as a refuge over the next century is inherently linked to projections about species' population status, because for some species the LIA will be the last remaining summer sea ice habitat e.g. ref. 27 . It is critical that future work quantify the resilience of this area for conservation and management of ice-dependent mammals under climate change.

Methods
Model and model configuration details. PIOMAS consists of coupled sea ice and ocean model components 28 . The sea ice model is a multi-category thickness and enthalpy distribution sea ice model which employs a teardrop viscous plastic rheology 29 , a mechanical redistribution function for ice ridging 30 PIOMAS is capable of assimilating satellite sea ice concentration data using an optimal interpolation approach 36 either over the whole ice-covered area or only near ice edge. In our run HIST, satellite ice concentrations are assimilated only near the ice edge (defined as 0.15 ice concentration). This means that no assimilation is conducted in the areas where both model and satellite ice concentrations are above 0.15. If the observed ice edge exceeds the model ice edge, then sea ice is added to the thinnest sea ice thickness category and sea surface temperature (SST) is set to the freezing point. If the model ice edge exceeds observations, excess ice is removed in all thickness categories proportionally. This ice-edge assimilation approach forces the simulated ice edge close to observations, while preventing satellitederived ice concentrations (which can be biased low during the summer e.g. ref. 37 ) from inaccurately correcting model ice concentrations in the interior of the ice  PIOMAS also assimilates SST 40 , using observations provided in the NCAR/ NCEP reanalysis (see below for atmospheric forcing) which in turn are derived from NOAA's OISSTv2.1 data set 41 . SST assimilation is only conducted in the open water areas, not in the ice-covered areas to avoid introducing an additional heat source into the sea ice budget 42 . For this study, we also conducted a number of sensitivity simulations in which no assimilation of ice concentration and SST is performed (see below).
Daily mean NCEP/NCAR reanalysis data 43 are used as atmospheric forcing, i.e., 10-m surface winds, 2-m surface air temperature, specific humidity, precipitation, evaporation, downwelling longwave radiation, sea level pressure, and cloud fraction. Cloud fraction is used to calculate downwelling shortwave radiation following Parkinson and Kellogg 44 . Precipitation less evaporation is calculated from precipitation and latent heat fluxes provided by the reanalysis model and specified at monthly time resolution to allow the calculation of snow depth over sea ice and input of fresh water into the ocean. There is no explicit representation of melt-ponds in this version of PIOMAS. River runoff into the model domain is specified from climatology 45 . Because of the uncertainty of net precipitation and river runoff, the surface ocean salinity is restored to a salinity climatology 46 with a 3-year restoring constant. Surface atmospheric momentum and turbulent heat fluxes are calculated using a surface layer model 47 that is part of the PIOMAS framework. Additional model information can be found in in Zhang and Rothrock 28 . PIOMAS has undergone substantial validation 23,[48][49][50][51][52] and has been shown to simulate sea ice thickness with error statistics similar to the uncertainty of the observations 52 . Validation results for ocean profiles for the WS are shown in S5.
Sea ice mass and upper ocean heat budgets. Components of the sea ice mass and upper ocean heat budgets are computed directly from model output and residuals. F prod is calculated as F prod = Δh/Δt -F adv and F bot = F prod -F atm-ice . All heat entering the uppermost ocean grid cell is used to melt ice until SIT = 0; however, subsurface shortwave radiation penetration and attenuation are allowed, which can warm the ocean below the uppermost grid cell. F ocndyn over the upper 60 m (Eq. (1)) can be partitioned into F ocndyn = F ocnadv + F diff + F convect where F ocnadv, F diff , and F convect are heat exchanges between the upper 60 m of the WS and the adjacent ocean via horizontal and vertical advection, horizonal and vertical diffusion, and vertical convection, respectively. F ocnadv is calculated directly from model ocean temperatures and velocities, and the sum of F diff + F convect found as a residual, i.e., F diff + F convect = ΔH/Δt -F atm-ocn − F bot − F ocnadv where ΔH/Δt is calculated directly from model temperature profiles. We find that F diff + F convect is negligible, meaning that horizontal and vertical advection terms (more formally, heat flux convergence) dominate. This is illustrated in Fig. S6, which shows a strong ocean warming within~100 km of the north Greenland coast owing to lateral heat flux convergence. This is nearly exactly balanced (not shown) by the vertical fluxes, i.e., downwelling, in keeping with previous results 53 . Finally, by comparing the heat budget for summer 2020 simulations with and without data assimilation (i.e., HIST vs. INIT), we find that this numerical effect produces only a negligible heat flux term and so is neglected here (it might be larger in other regions or over a longer time period of simulation). All ice mass and ocean heat budget terms are presented in units of meters of ice melt, assuming an ice density of 917 kg m −3 and latent heat of sea ice fusion of 3.293 × 10 5 J/kg.
HIST, INIT, and ATMOS Runs. The single HIST simulation uses data assimilation for the entire simulation period and is the basis for our analysis except for the sensitivity experiments described next. The INIT and ATMOS ensemble runs turn off the assimilation after May 31, 2020. For the JJA period of comparison, differences between the HIST run (which includes assimilation) and the equivalent members from the following ensemble runs (which do not include assimilation) are negligible.
Attribution of drivers. The INIT and ATMOS ensembles allow a partitioning of the proximate causes of the 2020 sea ice anomaly into those driven by the initial spring conditions (sea ice and ocean) and those related to the evolution of the atmosphere (winds, temperature, radiation, humidity) over the summer. To compute the relative contribution, we calculate spatially averaged SIC and SIT differences between the INIT and ATMOS ensemble medians and the HIST median at the time of the observed and simulated WS SIC sea ice minimum (August 14). The ensemble median here represents normal conditions as the reference to which conditions (sea ice for INIT, atmosphere for ATMOS) being tested are compared. The difference from HIST is considered the contribution of the respective 2020 condition, initial ice thickness for "INIT" and atmosphere for "ATMOS." This difference in SIC (SIT) is 6.3% (0.35 m) for INIT and 31% (1.5 m) for ATMOS. Adding these differences yields a total SIC response of 37.3%, and with respective fractions for "INIT and ATMOS" yields a 17% (6.3%/37.3%) role of initial conditions and a 83% (31%/37.3%) role for the atmosphere. The impact on SIT is slightly higher with respective contributions of 19% and 81%.
This partitioning can be used as a measure of the relative impacts of climate change and internal variability. Loosely following the framework of Trenberth at al. we assume atmospheric variability is governed by internal variability, and initial (i.e., spring) sea ice conditions to be driven by long-term climate change. Therefore the 20%/80% partitioning provides an approximate measure of the contributions of climate change and internal variability on the 2020 event. This separation is not perfect because atmospheric warming appears to be playing a role as evident in the fact that ATMOS ensemble members 2018/2019 both yield ice concentrations well below the 1979-2020 mean/median. The assumption that initial ice conditions are entirely due to climate change is also not entirely correct either, since internal variability also plays a role in sea ice conditions 16 . Nevertheless, our experiments clearly show that the climate signal of thinning sea ice exerts an impact on the magnitude of internally driven extreme events in the WS. Moreover, the fact that dynamic thickening of WS spring sea ice conditions (likely the result of internal variability) did little to improve the resilience of sea ice later in the summer provides an indication that climate change-driven thinning will likely influence future events.
Model uncertainties. As noted, PIOMAS has undergone substantial validation with respect to sea ice thickness, volume 23,48,50-52 and motion 54 . A measure of the uncertainty of ice-mass budget terms can be obtained from a recent study 55 that compared monthly advection and ice production terms from PIOMAS with another numerical model and estimates derived from satellite observations. Mass budget terms from the three different sources are highly correlated and provide confidence that the relationship of budget terms is correct even if their magnitudes may have error. In addition, our INIT and ATMOS model simulations incur additional uncertainties due to the lack of a direct feedback between the atmosphere and ice-ocean system. However, this problem is less severe in the summer season which is our focus here, because summer thermal contrasts are small between the marine surface and the atmosphere. Future experiments with coupled models that allow for a "replay" of observed variability will be needed to verify this.  56 . For 2020, we use the seaice_conc_cdr variable from the Near Real Time Version of the NSIDC-CDR 57 https://doi.org/10.7265/ N5FF3QJ6. This variable is derived using the identical algorithm as the "goddard_merged_seaice_conc" and therefore should provide a consistent record. Prior to 1987 the data are from SMMR and thus are only provided every other day and with a polar data gap (i.e., "pole hole") that straddles the northern edge of our domain. We, therefore, set ice concentrations to 100% in this area; the SIC minimum series is not affected by this filling.
Artist Sea Ice Concentrations (ASI). The ARTIST sea ice algorithm (ASI) 58 provides daily high-resolution ice concentration data derived from the 89 Ghz passive microwave channels on SSM/I and SSM/IS and from the AMSR-(E/2) sensors on board of Aqua and Global Change Observation Mission (GCOM). Data sets were downloaded from the University of Bremen, https://seaice.uni-bremen.de/data/amsr2/asi_daygrid_swath/ n3125. This data set has higher resolution relative to the NSIDC-CDR SIC's, but it only covers the period 2001-present. It thus complements the NSIDC-CDR SIC data set which offers the full 41-year record.
NSIDC Ice Age. This product provides the age of sea ice up to 15-16 years old. The age of sea ice is computed by tracking the motion of contiguous areas of sea ice from observed ice motion data derived by blending sea ice motion derived from passive microwave, visible band AVHRR data and in-situ observations 59 . https://doi.org/ 10.5067/UTAV7490FEPB.
AWI Cryos/SMOS. The Alfred Wegener Institute provides weekly and monthly products of sea ice thickness derived from the ESA CryoSat-2 radar altimeter merged with ice thickness derived from the SMOS L-band passive microwave instrument 60 . CS-2 provides a more accurate measurement for thicker ice, while SMOS provides ice thickness data for sea ice up to 1 m thickness. Neither product is available from May through September when increases in snow water content make retrievals of thickness too uncertain to be useful. Data were downloaded from ftp.awi.de:/sea_ice/product/ cryosat2_smos/v203/nh.  Figure S5 are from the Word Ocean Data Base https://www.ncei.noaa.gov/products/world-ocean-database 61 .