Rise and fall of ice production in the Arctic Ocean’s ice factories

,


Introduction
Arctic sea ice extent and thickness is in an on-going decline, sustained over at least the length of the satellite record, and directly related to anthropogenic carbon emissions. 1The retreat of Arctic sea ice is both an expression of and a driver of Arctic Amplification, which is the phenomenon of intensified climate change in the Arctic region relative to lower latitudes. 2,3,4 Tgh losses are observed in all seasons, they are more pronounced during late summer than late winter, in terms of both extent 5 and thickness. 6These trends reveal that increases in summer melting are to some extent compensated by increasing ice production during winter, 6 though not by enough to prevent the continued decline of the annual mean reservoir of sea ice.
As progressively more ice melts during the summer, winter ice production is critical in restoring the ice pack before the onset of polar day, when the high albedo of sea ice and its snow cover plays a crucial role in the radiative budget of the region.
The distribution of sea ice is also important during the winter; it limits heat fluxes from the ocean to atmosphere and plays a complex part in the surface momentum balance. 7Sea ice growth also plays an important hydrographic role: brine is rejected from sea ice during freezing, and freshwater is redistributed via sea ice motion, which acts as a vector between locations of growth and melt. 8,9 this paper, we aim to further understanding of the processes that control winter sea ice production in the Arctic, and explore the competing dynamics of these processes under climate change.
There is an apparent tension in the observation of rising winter sea ice growth in concert with Arctic Amplification, which is intensified during the winter. 3,10 owever, the two phenomena are most likely linked. 11During the autumn and winter, heat sequestered in the Arctic Ocean during the summer is released into the atmosphere prior to and during sea ice growth.Sea ice thinner than ∼0.4 m thick permits heat fluxes one to two orders of magnitude larger than those through perennial ice; 12 thickness declines during winter are an important part of Arctic Amplification, to which Lang et al. 13 attribute a rise of ∼1°C in Arctic surface air temperature per decade.In principle, this surface warming may be partially offset by increased long-wave emission. 14However, due to the predominantly stable atmospheric stratification, surface warming is inefficiently transmitted to the top of the atmosphere, and therefore outgoing longwave emissions only very weakly compensate the warming (the lapse-rate feedback). 15,16,17 Sice growth is inextricably linked to heat fluxes from the ocean to the atmosphere in the Arctic, because growth rate is determined by the energy balance at the lower boundary of the ice, the interface at which ice grows.As ice thickens, the conductive heat loss to the atmosphere declines, and the growth rate slows in direct proportion.Thin ice therefore not only permits higher heat fluxes across it, but also grows more rapidly than thicker ice.This thickness-growth relationship underpins a well-known negative feedback on Arctic sea ice loss. 18,19 here are other effects, however, that also act to promote sea ice growth in a warming Arctic.For one thing, ice that is formed later in the season is more likely to avoid the limiting effects of snow-an effective insulator-on growth rate. 20,21 imilarly, snow melt during summer preconditions more rapid growth at the beginning of the freezing season due to the loss of this insulating snow cover. 22Additionally, thinner sea ice is weaker and more mobile, leading to an increase in wind-driven sea ice divergence, 23,24,25 which has long been recognised to promote winter growth. 26These negative feedbacks provide stability to the Arctic sea ice system, such that rapid or irreversible losses for summer sea ice area are unlikely. 27,21,14 Iecent warm winters such as 2015/2016 and 2016/2017, growth suppression and even winter melting have been observed, 28,29,19,30 calling into question whether the action of these negative feedbacks may be becoming overwhelmed by warming.Based on simulations with the sea ice model CICE, Stroeve et al. 19 suggest that the thermodynamic growth during the winter 2016/2017 was 11-13 cm down on the 2011-2017 mean, with the overall positive trend in winter ice growth from 1985 ending in 2012; ice growth then weakens until the end of the model run in 2018.Using the same coupled climate model we employ here (CESM-LE), Petty et al. 31 identify that the positive correlation between temperature and winter growth weakens through the mid-century under the RCP8.5 emissions scenario, eventually becoming a negative relationship-indicative of warming overwhelming the negative feedbacks on sea ice loss.
Observational evidence of increasing winter sea ice production, and model-based suggestions of an imminent or recentlybegun decline-both associated with anthropogenic climate change-motivate the questions: when will this transition between rising and falling ice production begin?What physical indicators may help us to predict when we are at this turning point?
Building on the recent work of Petty et al. and Stroeve et al., can we understand the changes in ice production in terms of underlying physical processes?
We seek here to answer these questions with a regional focus on the Kara and Laptev seas, a region which has been referred to as the 'ice factory' of the Arctic Ocean 32 (Fig. 1a).Topologically, the Kara and Laptev seas are classified as interior shelves of the Arctic, 33 and are relatively fresh, receiving 50% of the freshwater runoff to the Arctic Ocean. 34During winter, winds predominantly drive sea ice northwards from the Siberian coastline in the Kara and Laptev seas (Fig. 1b), opening up perennial flaw leads/polynyas, where open water separates landfast ice from mobile pack ice.Sea ice forms rapidly in these flaw leads, and a number of studies have sought to constrain the outsized contribution to the Arctic sea ice budget made by these regions. 35,36,37,38,39,40 The ccurrence of wind patterns that increase the advection of sea ice away from the Siberian coasts thus intensifies sea ice production in these regions, and increases the fraction of (relatively thin) sea ice in the Arctic basin originating in the Kara and Laptev seas. 41,42 igh rates of freezing in these shallow seas are also thought to help bolster the cold halocline of the central Arctic basins. 43,44 he stability of the Arctic halocline has emerged as a key climate change indicator, 45 in light of episodic collapses of the winter halocline in the Eurasian basin and northern Barents Sea, 46, 47 which lead to the shoaling of Atlantic Waters and increased heat fluxes to the surface.We develop a simple linear model for ice production that is informed by the physics of sea ice growth, and can successfully explain both the rise and fall of ice production in this region under climate change.We compare the key climate variables from CESM-LE to equivalent variables from a range of observation-based products, and use these variables in our linear model to reconstruct ice production based on this linear model.The reconstruction suggests that, as in CESM-LE, we may presently be passing peak ice production.

Shrinking growth season, rising growth rate
During the summer melt season, large areas of open water develop in the Kara and Laptev seas as the sea ice edge retreats towards the pole (Fig. 2a).The winter refreeze, however, begins along the Siberian coastline. 50In CESM-LE, the newly formed ice along the Siberian coastline reconnects with the sea ice of the interior Arctic Ocean in an hourglass shape that becomes progressively pinched, and occurs later, through the 21st Century (Fig. 2).
On average, this refreezing in the early winter is the most productive period for winter ice growth (Fig. 3), and is linked to the areal expansion of sea ice (Fig. 2).After an initial peak, ice production in the ensemble mean gradually curtails through the winter-as ice and overlying snow cover thicken-before diminishing more steeply in April.As the climate warms in CESM-LE, however, October ice production rapidly falls, approaching zero by the 2030s (Fig. 3).(Note that ice production as plotted here only considers freezing, and not the negative contribution from melting.)Initially this trend is compensated by a rapid recovery and a higher sea ice production peak in November, connected to the freeze-back of a larger open water area (c.f.Fig. 2).The sea ice production peak itself shifts later after the 2010s, continuing to increase in magnitude until the 2040s, before falling rapidly to mid 20th Century levels by 2080.
The overall trend is one of progressively delayed freezing onset, followed by ice production that is more intense in the midwinter throughout RCP8.5 than it is in the 20th Century.The graph-area under the ice production curves in Figure 3 represents the total winter ice production; our quantity of interest.In the ensemble mean, ice production rises gently from approximately 1970 to 2010, before falling from 2020 onwards (Fig. 4a).
In examining the causes behind this rise and fall, we start by recognising that sea ice production is limited by the (locationdependent) duration of the freezing season, and the rate of growth where and when it is freezing.
Indeed, we can decompose total ice production into two components (Fig. 4).Firstly, the spatio-temporal duration of freezing, which can be expressed as the number of freezing area days (Fig. 4b).To extend the ice factory metaphor, this can be viewed as the usage of the ice factory: the area and time over which it is "switched on".Secondly, we can consider the mean growth rate over the regions where, and times when, sea ice is growing (Fig. 4c).In the ice factory metaphor, this is the efficiency of the ice factory when and where it is in use.We calculate the mean winter growth rate simply as the total winter ice production divided by the number of freezing area days.
The decomposition reveals quite different trends.The mean number of freezing area days (or ice factory usage), begins to decline from the 1990s onwards, which accelerates after 2050 (Fig. 4b).The mean growth rate (or ice factory efficiency), on In the absence of incoming shortwave radiation, the growth rate of ice is effectively captured by the balance of heat fluxes at its base. 51The growth of sea ice, of thickness h i , depends on a greater conductive heat loss to the atmosphere through the ice, F c than heat flux from the ocean at its base, F w , as shown in Equation 1, where ρ i is the density of ice, L i is the latent heat release/uptake on freezing/melting, and fluxes are positive upwards.
As sea ice thickens, the conductive heat fluxes F c through the ice diminish, bringing the growth rate down towards zero, depending on the size of the oceanic heat flux, F w .Other factors besides ice thickness also control F c , including the air temperature and snow cover.Assuming a layer of snow with thickness, h s overlying ice with thickness, h i , with corresponding thermal conductivities λ s and λ i , respectively, and linear temperature profiles through both layers, the conductive heat flux is 125 where k is an effective heat transfer coefficient between the snow/ice surface and the atmosphere, and ∆T = T w − SAT; in which T w is the freezing temperature of water and SAT is the surface atmospheric temperature.
In the Arctic Ocean, F w is generally low during the freezing season, of order 1 Wm −2 or less in the basin interiors, but sometimes higher directly over the Atlantic Water boundary current or over rough topography. 52Unlike in the Southern Ocean, 20 1960 2000

A linear model for ice production
By revealing that the number of freezing area days and mean winter growth rate exhibit quite different characteristics, the decomposition in Fig. 4 offers a high-level insight into the competing processes governing ice production changes.We now 135 seek to isolate a set of observable variables that can capture the physics of these processes.We then explain the changes in ice production through time with a linear model involving these variables.
As discussed above, variability in the number of freezing area days during Oct-Apr is dominantly captured by variation in the Sep sea surface temperature.The mean growth rate, meanwhile, must be controlled by the variables affecting the heat fluxes F w and F c (Equation 1).We can diagnose F w directly from the climate model and observations (though the latter are sparse), but we seek to capture F c in terms of readily observable climate variables, using insight from Equation 2.
The conductive heat flux is inversely dependent on snow and ice thicknesses (Equation 2).Since snow is much more insulating than ice (λ s ≈ 0.1λ i ), even a thin layer of snow may significantly affect the conductive heat flux. 51Equation 2 also shows that as surface air temperatures rise-and ∆T lessens-the conductive heat flux decreases in linear relation, as does the growth rate (Equation 1).Positive correlations between surface air temperatures and growth rate 31 rest on a causality in the opposite direction: increased growth rates and attendant upwards heat fluxes serve to boost surface air temperatures.We must therefore include the effect of air temperature in a way that explicitly captures its causal impact on sea ice growth.We do this-as detailed below-by mimicking the growth rate equation, and including ∆T as a factor in each of our regressors.
The sea ice thickness h i features on both sides of the growth rate equation (Equation 1) due to its involvement in F c .
Therefore, to avoid circular logic we must not include the evolution of h i in our linear model.However, we can consider initial conditions from which h i evolves according to the given growth rate.These initial conditions can be both the ice thickness at the start of the season, and occasions during the season when h i is set to zero (by divergence) as an initial condition for growth.The minimum sea ice area (usually occurring in September) determines the area over which h i = 0 at the start of the Divergence of the Arctic ice pack is spatially and temporally complex, and to understand its effect on thermodynamic growth we must also consider the role of convergence, which thickens ice through ridging and rafting, and may partially balance divergence that occurs in a given area.Because the growth rate of ice is an inverse function of its thickness (Equation 2), and because ice cracks under tension-thus setting h i = 0 in the diverged area-we can expect the effect of convergence and divergence on growth rate to be asymmetric.Indeed, in grid cells that exhibit both sea ice divergence and freezing within the Kara-Laptev region, divergence explains 37% of the variability in thermodynamic ice growth.In convergent-freezing cells on the other hand, there is no significant relationship.Because the effect of divergence on growth rate is not balanced by the effect of convergence, we must consider not only net divergence over the region, but also divergence that is balanced by convergence elsewhere within the region: we term this the compensated divergence.There is, by design, no double-counting between the net and compensated divergence.The two are therefore only weakly correlated-preferable for the purposes of building a linear model.
We We can include the effect of snow as ∆T /h s , as per Equation 2, where h s is the mean winter snow depth wherever ice is present.Gathering together these terms, our linear model for total winter ice production, including regression coefficients, β n is then: Where h s is the mean snow thickness where ice is present, A Sep is the open water area available in September, A net is the net area diverged integrated over the winter season, A comp is the compensated area diverged integrated over the winter season, SST Sep is the September sea surface temperature (the top 10 m ocean temperature), and ∆T is the temperature difference across the ice, determined by the surface atmospheric temperature only as we set T w to −1.8 • C (the freezing temperature).
We do not include the ocean-to-ice heat flux term, F w ; this is simply because, when included in this multiple regression, it is statistically insignificant.
We perform multiple linear regression across the ensemble members of CESM-LE to determine β n , using internal variability during a) the 20th Century run, b) the RCP8.5 run and c) the full timeseries (Fig. 5).We plot standardised coefficients, using 20th Century run standard deviation values.This indicates which regressors have the largest effect on ice production in terms of the natural range of variability in the 20th Century run.We preserve this standardisation in the plots for the RCP8.5 run and the full timeseries coefficients in order to compare absolute magnitudes.The variance in the CESM-LE internal variability explained by the linear model is consistently high: 81% during the 20th Century run, 76% during the RCP8.5 run and 78% over the full run.
There is relatively good agreement between the coefficients determined using the 20th Century, RCP8.5 and full timeseries of CESM-LE, and unanimous agreement in the sign.This suggests that the linear model is robust to different climate conditions.The absolute magnitude of the ∆T /h s coefficient (β 1 )is remarkably similar between runs (Fig. 5).We can estimate the uncertainty associated with these estimates by limiting the dataset to individual ensemble members (grey points, Fig. 5).In the main, individual ensemble members yield estimates that cluster around the full ensemble value.However, for the coefficient associated with September open water ×∆T in RCP8.5, the full ensemble value lies outside this range; an effect likely due to September open water area becoming an increasingly skewed variable with time, as it is limited to the total area under consideration.
Three regressors are consistently dominant in their contribution to the internal variability: net divergence ×∆T , Sep open water ×∆T , and inverse snow depth ×∆T .It is worth noting, however, that their role in forced changes in ice production-as given by the ensemble mean-also depends on the magnitude of forced trends in each regressor.

Explaining forced changes in ice production
Having seen that our linear model can explain a large fraction of the internal variability in sea ice production in CESM-LE, we now ask whether it can explain the forced changes in sea ice production associated with climate change (Fig. 6a) and, in Bars show the regression coefficients calculated using internal variability in all 35 ensemble members, with the ensemble mean removed.
Variance explained is as follows.20C, R 2 = 0.81; RCP8.5, R 2 = 0.76; Full, R 2 = 0.78.Uncertainty is indicated by the spread of grey points, which show the regression coefficients calculated using each individual ensemble member.
doing so, help us to understand the process at play.To do this, we construct a timeseries of ice production using the regression coefficients calculated from internal variability in a) 20th Century data only; b) RCP8.5 data only and c) the full timeseries, combined with timeseries of the regressor climate variables from the ensemble mean.All three solutions for sea ice production successfully reconstruct the overall shape and interannual variability in the ensemble mean over the full period.The linear model also accurately captures the timing of the turning point from rising to falling ice production.
In light of the success of the linear model at reconstructing forced changes in sea ice production, we can next interrogate the contributions made by the different explanatory variables (Fig. 6b).Interpretation of these trends should critically consider that each regressor, except for Sep 10 m T, contains the product of surface air temperatures (as ∆T ) and another climate variable.
The trends of the underlying timeseries (Fig. 7) provide useful context.
The greatest contribution to increasing ice production up to 2020 comes from the β 2 ∆T A Sep term.This term can be interpreted as ice produced in the open water available at the start of the freezing season.As the September open water area increases during 1970-2020 (Figs.2a & 7b), this term contributes increasing ice production, despite SAT warming.As the ensemble mean Sep sea ice area approaches zero at around 2020 (Fig. 7b), this negative feedback on sea ice loss approaches Decreasing snow thickness on sea ice also provides an important contribution to rising ice production between c. 1980-2010 (Fig. 6b).While snow continues to thin to 2080 (Fig. 7f), after c. 2030 this thinning is out-competed by the effect of warming SAT, and the β 1 ∆T /h s term contributes to a decline in ice production.

220
The two terms associated with divergence show more muted forced changes, owing to a close balance through time between the positive effects of divergence and the negative effects of SAT on ice production.However, in both terms, warming air temperatures begin to out-compete increasing divergence from the mid-21st century onwards.The β 3 ∆T A net term also clearly contributes to ensemble mean interannual variability in ice production, reinforcing the important role of variability in winddriven sea ice divergence in explaining year-to-year changes in ice production in the Kara-Laptev region.The β 5 SST Sep term, though showing the weakest regression coefficient in the internal variability space, is the dominant driver of the forced decline in ice production.Indeed, Sep 10 m temperature exhibits a clear warming trend in the RCP8.5 run that rapidly departs from the range of 20th Century internal variability (Fig. 7e).There is evidently, then, a strong role for ocean heat derived from summer heating in delaying the freezing season and thus decreasing ice production.But is there a role for ocean heat in changing growth rate or ice production once the freezing season has begun?Ocean-to-ice heat fluxes (F w ) in the region are generally low during 20C at c. 0.5 W/m 2 , rising throughout RCP8.5 to nearly 2 W/m 2 by 2080.This, however, does not appear to noticeably impede ice growth in CESM-LE: F w is a statistically insignificant term when included in linear models for both ice production and growth rate, even in RCP8.5 (not shown).

Application to the observed Arctic
We next attempt to reconstruct changes in the Kara-Laptev region by applying observation-based estimates of the relevant climate variables to our linear model for ice production (Fig. 8).We use the regression coefficients derived from internal variability over the full run of CESM-LE.Note that this reconstruction using historical climate data is more analogous to the ice production in a single ensemble member of CESM-LE, rather than the (smooth) ensemble mean.The 10-year running mean in Unlike in CESM-LE, the β 2 ∆T A Sep term is not the dominant driver of the reconstructed increase in ice production.We can reconcile this by noticing the anticorrelation (R = −0.77) between the Sep SIA and SAT timeseries (Fig. 7a,b).Both timeseries show marked jumps in the winter 2005/2006-these changes act in opposing directions in terms of ice production and thus tend to cancel one another out.Given that the September open water area has already reached a maximum (as September SIA has declined to zero in recent years), we can expect future atmospheric warming to yield a negative ice production trend from this term, as per Fig. 6b.
Examining the trends from the other terms, we see that September SSTs exert a negative influence on ice production from the 2000s.The inverse snow thickness term β 1 ∆T /h s causes marked interannual variability and a strong positive trend in ice production: snow thinning out-competes warming over this period.
The contributions to ice production from the two divergent settings exhibit perhaps the strongest rises in the mid-part of the timeseries, out-competing atmospheric warming, until plateauing or shallowly declining in the final c. 15 years.

Discussion
Although CESM-LE is a complex coupled climate model, we have shown that we can understand variability and trends in ice production in the Arctic Ocean's "ice factories" using a simple equation.Our analysis shows that a number of negative feedback to ice loss-increasing September open water, increasing divergence, reducing snow depth-contribute to a gradual rise in ice production in CESM-LE from c. 1970-2010, but are increasingly outcompeted by atmospheric warming through the 21st Century under the RCP8.5 emissions pathway.Together with an important contribution from increasing upper ocean temperatures at the end of summer, this leads to a decline in ice production.

260
In CESM-LE, the timing of peak ice production is primarily set by the timing at which September sea ice area in the Kara-Laptev region approaches zero.At this point, one of the leading negative feedbacks to the loss of ice-the expansion of the September open water area-has reached a limit and can no longer contribute additional ice production.Continued atmospheric warming then ensures a decrease in the ice production attributed to this setting.In observational data, we see that the expansion of the September open water area has also reached a limit in recent years.As such, we should expect ice production in the Kara

265
and Laptev seas to decline in the coming decades under continued greenhouse forcing.
relevant climate variables, we see a similar trend to that shown by CESM-LE over the same time window.Reconstructed ice production appears to be currently passing a peak, with the recent warm winters of 2015/2016 and 2016/2017 28,29,19 notably triggering relatively low ice production.
Sea ice from the Kara and Laptev seas carries sediment, pollutants, trace elements and gases into the central Arctic and beyond. 54,35,55 Flling sea ice production in the region will therefore affect redistribution of biogeochemical matter in the Arctic Ocean, with implications for primary production and biodiversity. 54 the sea ice production in the Kara and Laptev seas begins to fall, we might question how other areas of the Arctic Ocean will fare.While much of the process-based understanding we have developed for the Kara-Laptev region is transferable to other areas, there are important differences in the mean sea ice state and oceanic setting that affect the balance of F c and F w .
As ice is generally thin in the Kara and Laptev seas, owing to the wind-driven time-mean divergence of sea ice, and because little ice survives the summer, F c is comparatively large (Equation 2).On the other hand, once the excess heat stored in the upper ocean has been removed to the atmosphere, F w is small.In other regions this balance may differ.
The same processes identified in our linear model that drive increasing ice production (by increasing F c ) in the Kara-Laptev region during 1970-2010 are likely to become increasingly important in the central Arctic.As the September sea ice edge progressively retreats through the central Arctic (Fig. 2a), we can expect year-on-year increases in Sep open water area and sea ice divergence, and decreases in snow depth in these regions-all of which drive increasing ice production and are negative feedbacks on sea ice loss.In CESM-LE, the share of Kara-Laptev ice production in the Arctic Ocean (excluding Barents Sea) declines from 25% in the 20th Century run to 19% by 2080 (it comprises 16% of the area).
These shifting patterns of ice production raise important questions for the stability of the Arctic halocline.High rates of freezing in the Kara and Laptev seas strengthen the cold halocline of the central Arctic basins, which limits upwards heat fluxes from Atlantic Waters and restricts the depth of water that is convectively cooled prior to and during sea ice growth. 56eezing in the Kara and Laptev seas can stabilise the halocline via the advective interleaving of cold shelf waters, densified by brine rejection, into the halocline of the Eurasian Basin. 43,57 nother mechanism is that winter convection down to the halocline is followed by surface injection of cold and fresh waters (from the melting of shelf-derived sea ice), thus renewing the cold halocline via a 'convective mode'. 44,58 oth mechanisms may operate in different seasonal and spatial contexts, 59,60,37 and both depend on ice production on the shelves.Meanwhile, increasing ice production and brine rejection in the Eurasian basin interior may weaken the halocline by encouraging vigorous convection, allowing heat from the Atlantic Waters to reach the surface. 48,46 ters deeper than the shelfbreak, we expect a more important role for F w in determining ice production on account of the availability of deep Atlantic heat.In the future, ocean-to-ice heat fluxes are likely to play a greater role in limiting sea ice production in the Arctic, a trend that may have already begun. 61In the shallow waters of the Kara and Laptev polynyas, on the other hand, winter convection may propagate to the seafloor, 32, 39 but because Atlantic Waters rarely upwell into these shallow areas where freezing is most intense, increases in F w from such mixing should be limited, though Atlantic Waters may warm halocline waters that intrude onto the shelf. 62Nonetheless, in our linear model for sea ice production, F w is an insignificant term even in the RCP8.5 run up to 2080.However, in other areas of the Arctic, and indeed the Antarctic, F w is likely to be important.
Our reconstruction of forced changes in sea ice production in the Kara-Laptev region (Fig. 6) shows that it is the heat stored in the upper ocean at the end of summer, rather than heat fluxed from below during the freezing season, that is more important in restricting ice production.A greater reservoir of summer heat requires more time to cool, thus shortening the freezing season (reducing the usage of the ice factory).September upper ocean temperatures can be expected to continue to rise as a result of increasing solar radiation (due to reducing summer sea ice concentration), increasing radiative forcing, and increasing (and warming) river runoff. 63,64 this paper we developed a simple linear model that captures ice production in key "ice factories" of the Arctic Ocean.
Ice production in these regions is most likely presently passing a peak, and will soon begin a long-term decline if greenhouse gas concentrations continue to rise.Reduced ice production in the Kara and Laptev seas will alter the physical and biological makeup of the Arctic Ocean, and we suggest that future work should address the profound impact of these forseeable changes.

Climate model
To high-emissions RCP8.5 run, 68 leading to over 4 • C warming by 2100.The RCP8.5 pathway at present slightly underestimates observed emissions; it presently represents a realistic or even conservative scenario. 69,70 dditionally, it provides big signal, which helps us to interpret climate dynamics in the presence of noise.The ensemble spread is entirely generated by simulated internal climate variability originating from very small, random differences in the initial air temperature fields.
Here we rely on monthly data between 1920-2080-sufficient time to assess major forced changes in ice production in the Kara and Laptev seas-from 35 ensemble members.Previous Arctic studies have investigated the impact of subsampling the ensemble. 71,72,73 Te number of ensemble members required varies according to the nature of the climate variables under investigation.For coarse metrics such as winter means, we are confident that 35 ensemble members is sufficient.
CESM1.1 and CESM-LE have been used in a broad swathe of Arctic climate studies and generally perform well. 74,71,75,76,77,73,78,79,31 Arctic  ice thickness in CESM-LE broadly corresponds in spatial mean pattern and trend during 1980-2015 with PIOMAS, 80 though exhibits slightly thicker (order several 10 cm) sea ice than PIOMAS.81,31 Although the satellite record of sea ice thickness is short, CESM-LE compares well with estimates of inter-seasonal thickness changes and interannual variability.31 Declining Arctic sea ice concentration and extent in CESM-LE has been more comprehensively compared to observations, which fit within its ensemble spread and compare well across all seasons.77,71,73 The changing open water season also shows correspondence with satellite observations, including in the Laptev Sea.73

Climate variables computed
To capture winter thermodynamic sea ice production, we use monthly data from October to April, inclusive, as per Petty et al. 31 In contrast with Petty et al., we analyse freshwater exchanges with the ocean, rather than changes in ice thickness, to explicitly isolate thermodynamic changes.Sea ice production is calculated as the sum of only the ocean-to-ice part of thermodynamic ice-ocean freshwater fluxes.
Our calculation of freezing area days using CESM-LE data is an approximation, dependent on the spatial discretisation of the ocean and sea ice components (nominally 1 • ×1 • ), temporal averaging (monthly) and number of ensemble members used (35).The mean winter growth rate in all such freezing grid cells is given by the total winter ice production divided by the number of freezing area days.
The other variables we consider here are: snow depth on sea ice (note that this is different to snow volume per unit grid cell area), sea ice concentration, sea ice divergence, surface air temperatures, and ocean temperatures in the top (10 m thick) grid cell.We investigate two types of sea ice divergence.Firstly, the net divergence, which-by the divergence theoremis equivalent to the area expansion of the sea ice area, or, when that area encompasses the whole Kara-Laptev region, it is the area export of sea ice from the region.We are also interested in divergence within the region that is compensated by convergence somewhere else in the region; this divergence can open or expand leads but still not lead to export.To capture this compensated divergence, we record the positive divergence, which is counted only over cells where the divergence is positive.
The compensated divergence is then found as compensated div = positive div − net div.
When plotting data, we identify the winter spanning years n and n + 1 with year n + 1; the winter 2015/2016 is labelled 2016.The variables extracted from the observation-derived estimates are computed in a consistent manner to their equivalents from CESM-LE.

Observation-based estimates
We use monthly mean 2 m temperatures from the ERA5 global reanalysis, 82 gridded at 0.25 • ×0.25 • resolution.ERA5 is the fifth generation of atmospheric reanalyses from the European Centre for Medium-Term Weather Forecasting (ECMWF).We use the terms surface air temperature and 2 m temperature interchangeably here.Surface air temperatures from ERA5 generally perform well in the Arctic relative to in situ observations, however consistently show warm biases of order several degrees over sea ice in winter months. 83Our analysis of these data spans the winters 1979/1980 to 2019/2020.
To plot mean sea ice velocities and calculate sea ice divergence, we make use of the National Snow and Ice Data Center (NSIDC) Polar Pathfinder v4.1 product, gridded at 25×25 km resolution and in weekly means. 49We analyse data from Polar Pathfinder encompassing the winters 1979/1980 to 2019/2020.
To document changes in sea ice area, we use sea ice concentration data from passive microwave satellite retrievals, prepared by the NASA Goddard Space Flight Centre (NASA GSFC).We use the 'merged' data, which combines data derived from the Bootstrap and Nasa Team algorithms, 84 and compares well with other estimates. 85,86 e use data spanning the winters 48

Figure 1 .
Figure 1.(a) Winter (Oct-Apr) mean ice production for the period 1979-2020 in the ensemble mean of CESM-LE.The Kara and Laptev seas study region is outlined in blue.(b) Mean sea ice velocities during 1979-2020 as recorded by the Polar Pathfinder v4.1 product, gridded at 25×25 km resolution and in weekly means. 49Means taken over periods when sea ice present.Shading indicates the magnitude of the mean velocities.Arrows are omitted in the Barents Sea, Baffin Bay and Greenland Sea, where the colourscale saturates, and flow is generally southwards.

Figure 2 .Figure 3 .
Figure 2. Decadal ensemble mean winter sea ice extent edges in CESM-LE from 1970-2080.Study region outlined in black.Ice edges defined at 15% concentration.Grey line marks the mean 2020s sea ice extent edge.

Figure 4 .
Figure 4. Winter ice production (a) and decomposition into freezing area days (b) and growth rate (c).Note that the freezing area days is rescaled in terms of the total area of the study region, such that the upper limit is 210 on the y axis: the duration of our seven month winter period in days.
season: the September open water area.In addition, open water is created near-continually throughout the freezing season by divergence.Under tension, ice cracks in a brittle manner rather than thinning; the open water area is equal to the area diverged.
thus classify three settings in which open water is made available for ice growth: the September open water, created by summer melting and ice advection; the open water area due to net divergence; and the open water area due to compensated divergence.The sea ice growth in such open water areas will be linearly related to ∆T and the area that has been opened.In our linear model, we neglect the timing of divergence and simply use the total area diverged during the winter season in both cases.

Figure 5 .
Figure 5. Winter ice production regression coefficients, standardised by the respective one standard deviation values in the 20th Century run.

Figure 6 .
Figure 6.a) Ensemble mean winter sea ice production (black line), and reconstruction of ensemble mean using regressions trained on internal variability over the 20th Century run (light blue line), RCP8.5 run (light orange line) and full 1920-2080 timeseries (heavy blue line).b) Contribution to the reconstruction made by each component of the linear model, when trained upon internal variability from the full timeseries, plotted relative to the first winter (1920/1921). 35

Figure 7 .
Figure 7. Timeseries of linear model inputs from CESM-LE, with real-world estimates superimposed in red: a) surface air temperatures, b) September sea ice area, c) net area diverged, d) compensated area diverged, e) September 10 m T in CESM-LE and SST in NOAA OI SSTv2, f) snow depth on ice.The 20th Century run is shown in blue and the RCP8.5 run in orange; ensemble mean indicated by darker line.

Fig. 8 ,
Fig. 8, however, exhibits a strikingly similarity to the ensemble mean trend in CESM-LE: reconstructed winter ice production increases by 150-200 km 3 from 1983 to 2000.From about the mid-2000s, ice production appears to show a slight decline.While uncertainty is difficult to constrain from the underlying observation-based estimates, we note that the reconstruction is somewhat sensitive to the choice of product for SST, especially from 2006/2007 onwards.

Figure 8 .
Figure 8. Reconstruction of historical changes in the winter ice production in the Kara and Laptev seas relative to the first winter of the reconstruction (1982/1983).Reconstruction uses climate variables from observation-based estimates and the linear model trained on CESM-LE internal variability.Data are plotted with respect to the new year of the winter in question.Dashed line is the result of using NOAA OI V2 SST product, while the solid line uses HadISST.
analyse the forced response of sea ice production in the Kara and Laptev seas to climate change, we use data from 35 of the 40 ensemble members of the Community Earth System Model Large Ensemble (CESM-LE). 65A large ensemble such as CESM-LE permits investigation of forced responses to climate change in the context of internal climate variability. 65CESM-LE is based on the fully coupled model CESM1.1, and comprises the Community Atmosphere Model, version 5 (CAM5); the Parallel Ocean Program, version 2 (POP2); the Community Land Model, version 4; and the Los Alamos Sea Ice Model CICE, version 4 as its sea ice component.The version of CICE used features improvements to shortwave radiation interactions, including the effects of melt ponds and aerosol deposition on ice. 66The spatial resolution of the CESM1.1 ocean and sea ice models is nominally 1 • ×1 • longitude by latitude, while the atmospheric model is 0.9 • ×1.25 • .All CESM-LE ensemble members are forced by the same external forcing data, separated into two runs.Firstly, a run from 1920-2005 forced by historical external forcing data 67 that we refer to as the 20th Century run; secondly, from 2006-2100, a 1978/1979 to 2018/2019.For snow depth, we use data from a Lagrangian snow evolution model, SnowModel-LG, which provides daily data on a 25 km × 25 km grid and is forced by the MERRA-2 atmospheric reanalysis.87SnowModel-LG compares well to observational datasets in spatial and seasonal variability of snow depth and density.88The dataset spans winters 1981/1982 to 2018/2019.Consistent with our snow depth acquisitions from CESM-LE, we present snow thickness on sea ice where present.For sea surface temperatures, we use the Met Office Hadley Centre's sea ice and sea surface temperature (SST) data set, HadISST1.89The monthly data are provided on a 1 • ×1 • grid.We use data spanning the winters 1961/1962 to 2019/2020.The SST data are taken from the Met Office Marine Data Bank (MDB), as well as data from the Comprehensive Ocean-Atmosphere Data Set (COADS) (now ICOADS) where there were no MDB data.The sea ice data are taken from a variety of sources including passive microwave retrievals and digitized sea ice charts.Additionally, we use monthly-mean SST data at a 1 • ×1 • resolution from the National Oceanic and Atmospheric Administration Optimum Interpolation Sea Surface Temperature v2 (NOAA OISSTv2) dataset.90We use NOAA OISSTv2 data from the winters 1981/1982 to 2019/2020.