Internal Ocean Dynamics Control the Long-Term Evolution of Weddell Sea Polynya Activity

Open-ocean polynyas effectively couple the ocean and atmosphere through large ice-free areas within the sea-ice cover, release vast quantities of oceanic heat, and impact deep ocean ventilation. Changes in polynya activity, particularly in the Weddell Sea, may be key to longer time-scale climate fluctuations, feedbacks and abrupt change. While changes in the occurrence of Weddell Sea polynyas are generally attributed to changes in the atmospheric surface forcing, the role of internal ocean dynamics for polynya variability is not well-resolved. In this study we employ a global coupled ocean-sea ice model with a repeating annual atmospheric cycle to explore changes in Weddell Sea water mass properties, stratification and ocean circulation driven by open-ocean polynyas. During the 1300-year long simulation, two large polynyas occur in the central Weddell Sea. Our results suggest that Weddell polynyas may be triggered without inter-annual changes in the atmospheric forcing. This highlights the role of ocean processes in preconditioning and triggering open-ocean polynyas on multi-centennial time-scales. The simulated polynyas form due to internal ocean-sea ice dynamics associated with a slow build-up and subsequent release of subsurface heat. A strong stratification and weak vertical mixing is necessary for building the subsurface heat reservoir. Once the water column turns unstable, enhanced vertical mixing of warm and saline waters into the surface layer causes efficient sea ice melt and the polynya appears. Subsequent, vigorous deep convection is maintained through upwelling of warm deep water leading to enhanced bottom water formation. We find a cessation of simulated deep convection and polynya activity due to long-term cooling and freshening of the subsurface heat reservoir. As subsurface waters in the Southern Ocean are now becoming warmer and saltier, we speculate that larger and more persistent Weddell polynyas could become more frequent in the future.


INTRODUCTION
The Southern Ocean is an important region for formation of cold and dense bottom waters feeding the global ocean circulation and plays a key role in regulating Earth's climate by redistributing heat and regulating the uptake of atmospheric CO 2 (Orsi et al., 1999;Marshall and Speer, 2012;Watson et al., 2014). In the Weddell Sea, bottom water is formed on the shallow continental shelf in small coastal polynyas, as sea ice is pushed away from the Antarctic continent by strong winds, exposing the relatively warm ocean to the cold polar atmosphere. Intense sea ice formation and brine rejection leads to the formation of high salinity shelf water that spills over and forms Antarctic Bottom Water (AABW) (Orsi et al., 1999). This process, however, is poorly represented in today's climate models with the majority of the models contributing to the 5th phase of the Coupled Model Intercomparison Project (CMIP5) forming bottom water by open-ocean deep convection in the Weddell Sea (Heuze et al., 2013;Heywood et al., 2014). This often leads to significant biases in modeled Southern Ocean stratification affecting deep ocean volume transports (Timmermann and Beckmann, 2004;Behrens, 2016). This highlights the importance of understanding the processes leading to deep convection in climate models in order to improve future projections (Zhang et al., 2019).
In the models, deep convection is closely linked to the formation of large open-ocean (or sensible heat) polynyas in the Weddell Sea. As opposed to coastal polynyas, Weddell polynyas form off-shore in response to upwelling of Warm Deep Water (WDW) located below a weakly stratified surface layer (e.g., Martin et al., 2013;Cheon et al., 2015). In the winter, when the stratification is weak, small changes in surface buoyancy can trigger deep convection causing WDW to be mixed into the surface layer where it melts sea ice, forming a region of open water within the sea ice pack. Intense cooling at the air-sea interface forces surface water to cool and sink in the polynya, presenting an alternative source of deep and bottom waters (Killworth, 1983). Consequently, Weddell polynyas have widespread implications for deep ocean ventilation and uptake of anthropogenic and natural CO 2 Menviel et al., 2018), overturning circulation strength (Martin et al., 2013(Martin et al., , 2015, and impact Antarctic as well as global climate variability (Pedro et al., 2016;Cabré et al., 2017;Zhang et al., 2019). Recently, Naughten et al. (2019) also showed that Weddell Sea polynyas may accelerate melting of nearby Antarctic ice shelves by increasing the transport of WDW into the ice shelf cavity. Polynyas have also been invoked as a potential mechanism to explain rapid climate fluctuations in the North Atlantic during the last glacial period (Vettoretti and Peltier, 2016).
In contrast to the CMIP5 models, exhibiting frequent openocean deep convection, the occurrence of large polynyas is a rare phenomenon in the instrumental record. During the mid-1970s, satellite observations revealed a large (250,000-300,000 km 2 ) winter-persistent polynya in the eastern Weddell Sea (Carsey, 1980). This feature, which is known as the Great Weddell Polynya, persisted for three consecutive winters (1974)(1975)(1976) and formed after a prolonged period of cold and dry atmospheric conditions. This contributed to a saltier surface layer and a weakening of the surface stratification enabling WDW to reach the sea ice (Gordon et al., 2007). Several smaller and intermittent polynyas were observed during the 1990s over Maud Rise (Smedsrud, 2005), and again in 2016 and 2017 reaching a size of 50,000 km 2 (Swart et al., 2018;Campbell et al., 2019;Jena et al., 2019).
Unfortunately, due to the short instrumental record, there are large uncertainties in the long-term evolution of Southern Ocean polynya activity. Although recent attempts have been made to reconstruct past polynya activity, the results are still ambiguous (Goosse et al., 2021), and the climatic conditions necessary to sustain deep convection in the Southern Ocean are not well constrained. de Lavergne et al. (2014) speculated that open-ocean convection in polynyas may have presented a dominant mode of deep water formation in the past, but has reduced over the past couple of decades due to a freshening of the Southern Ocean surface waters. Similarly, other studies point to large-scale changes in atmospheric circulation patterns as being the controlling factor for polynya formation and its absence in recent decades (e.g., Gordon et al., 2007;Cheon et al., 2017;Campbell et al., 2019;Kaufman et al., 2020). Meanwhile, less attention has been given to how internal ocean-sea ice processes might influence stratification (especially over long timescales) and its role in preconditioning and potentially triggering openocean polynyas. In particular, given its importance for polynya formation, changes in WDW heat content (e.g., by mixing processes in the ocean) may present a major control on the frequency of polynya events. In climate models, the simulated open-ocean convection is often preceded by a build-up of heat and salt at mid-depths over several years or decades (Martin et al., 2013;Cheon et al., 2015;Zanowski et al., 2015;Dufour et al., 2017). The build-up eventually destabilizes the water column, triggering cycles of deep convection and subsequent build-up. This is also evident in observations by Robertson et al. (2002) and Smedsrud (2005), showing that the WDW warmed from the late 1970s to the 1990s when the Weddell Polynya was absent. Moreover, it has speculated that the observed warming of Southern Ocean deep waters and shrinking of AABW in the last decades may reflect changes in polynya activity (e.g., Purkey et al., 2012;Zanowski et al., 2015;Wang et al., 2016). While it remains debated whether subsurface heat accumulation, which is simulated in models, can trigger polynyas (e.g., Campbell et al., 2019), it nevertheless suggests a close connection between the WDW heat content and the occurrence of Weddell polynyas.
In this study, we explore the impact of long-term changes in WDW and Weddell Sea stratification and how it affects the formation of open-ocean polynyas. Applying a stand-alone ocean-sea ice model with a constant atmospheric forcing, we analyse a 1300-year transient simulation where two large Weddell polynyas form. This presents a unique opportunity to test how interior ocean processes influence the frequency of Southern Ocean deep convection events in the absence of atmospheric variability. Rather than trying to explain the trigger mechanism for the observed Weddell Polynya in the 1970s, this study focuses on the general dynamics of polynya formation by addressing two main questions: (1) Can build-up of subsurface heat and salt alone trigger Weddell polynyas? and (2) What controls the timescale and frequency of polynya events in the Southern Ocean?

Model Description
To study the polynya dynamics, we use a stand-alone oceansea ice configuration of the Norwegian Earth System Model (NorESM-OC1.2; Schwinger et al., 2016). The ocean-sea ice component of this NorESM-OC1.2 is the same as in Bentsen et al. (2013), and there is also a carbon-cycle part which has not been used in this study. Changes in the carbon-cycle induced by Southern Ocean convection is beyond the scope of this paper, but could be relevant to explore in future research.
The family of NorESM models is based on the Community Climate System Model version 4 (CCSM4), using the same land and sea-ice models, but with a different atmosphere and ocean component. The ocean model is based on the Miami Isopycnal Coordinate Ocean Model (MICOM). MICOM uses potential density as the vertical coordinate which differs significantly from z-coordinate models by mimicking the real structure of the ocean as discrete layers of constant density. One of the key advantages of the isopycnal coordinate models is reduced spurious mixing due to the numerical representation of isopycnal advection and mixing. This allows for a more accurate transport, improved representation of dense water overflows in gravity currents and preservation of water mass characteristics (Griffies et al., 2000;Assmann et al., 2010;Bentsen et al., 2013). The sea ice model is the Community Ice Code version 4 (CICE4) which is the same sea-ice component used in CCSM4 (Hunke and Lipscomb, 2008;Gent et al., 2011). CICE4 is a fully dynamicthermodynamic sea ice model and shares the same horizontal grid as the ocean component.
The model is configured on a tripolar grid with a longitudinal resolution of 2 • and a variable latitudinal grid spacing from 0.5 • at the Equator to 0.35 • at high southern latitudes. In the vertical, the ocean is divided into 51 isopycnal layers with potential density classes ranging from 23.602 to 33.2 kg m −3 . The first two model layers constitute the mixed layer where the density can evolve freely. The uppermost layer is limited to 10 m if the total mixed layer depth is <20 m, thus allowing the surface ocean to respond faster to surface fluxes. The mixed layer depth is parameterized according to the turbulent kinetic energy (TKE) model based on Oberhuber (1993). To reduce SSS biases in high latitude regions, salt released due to sea ice formation is evenly distributed below the mixed layer down to a depth where the density difference is 0.4 kg m −3 relative to the surface. This avoids build-up of salt and improves stratification in seasonally ice-covered oceans (Bentsen et al., 2013).
In isopycnal models potential density is referenced to a given pressure. This means that the flow characteristics are most accurately represented near the reference level, where isopycnal and neutral density surfaces are similar. When the pressure deviates significantly from the reference pressure, i.e., in the deep ocean, the potential density becomes increasingly nonneutral. In order to avoid this non-neutrality, most isopycnal models now use a potential density referenced to 2,000 db, which maximizes the neutrality of isopycnal layers (McDougall et al., 2005). In high-latitude regions, where the water column is weakly stratified, the stratification is more sensitive to the choice of the reference pressure. In a previous version of NorESM, Assmann et al. (2010) showed that the 2,000 db reference pressure caused the stratification in the Weddell Sea to be unstable, leading to overestimated vertical mixing which eroded the WDW. As the subsurface heat reservoir is essential for the dynamics of polynya formation and deep convection in the Weddell Sea (Martin et al., 2013;Gordon, 2014;Cheon et al., 2015;Dufour et al., 2017) we chose a reference pressure of 1,000 db. This ensures a stable stratification in the Weddell Sea region and improves the representation of the WDW layer. We note, however, that this could introduce biases in the deepest part of the basin where the pressure can differ significantly from the potential density at 1,000 db.
The isopycnal formulation used offers complete control of the diapycnal mixing applied in the model. The total diffusivity of diapycnal mixing consists of four components relating to different mixing processes in the interior ocean: (2) where κ b is a constant background diapycnal diffusivity set to 10 −5 m 2 s −1 , including a latitude dependency according to Gregg et al. (2003). κ r is the shear driven diapycnal mixing as a function of the local gradient Richardson number which relates to the balance between shear induced turbulence and buoyancy forcing. The tidally driven mixing κ t induced by the internal wave field and braking of internal waves is implemented from Simmons et al. (2004). Lastly, κ N relates to the diapycnal mixing when local stability is weak and is triggered if the local interface density difference is <0.006 kg m −3 . A convective adjustment is then made at each time step by increasing the diapycnal diffusivity to 0.05 m 2 s −1 . Mixing along isopycnal surfaces is parameterized according to Eden and Greatbatch (2008).

Experimental Setup
The model is forced with the Coordinated Ocean-ice Reference Experiment version 1 (CORE-I) repeated normal year forcing (CORE-NYF) which has been derived from 43 years of interannually varying forcing with a seamless transition from 31 December to 1 January (Large and Yeager, 2004). The CORE framework is used as a tool to study the behavior in coupled global ocean and sea ice models forced by a common atmospheric dataset (Griffies et al., 2009). This reduces the complexity of the coupled climate system by prescribing atmospheric boundary conditions without dealing with the ocean-atmosphere coupling and feedbacks. The normal year forcing is constructed by repeating a 1-year annual cycle, which means that there is no variability in the atmospheric forcing beyond seasonal variations. Thus, the inter-annual variability simulated in the coupled oceansea ice model is due to internal dynamics. This approach allows us to explore the physical processes internal to the ocean-sea ice system in response to a repeated annual atmospheric cycle and test if Weddell Sea polynyas can form without invoking changes to the atmospheric forcing. The model is initialized with January mean temperature and salinity fields from the Polar Science Center Hydrographic Climatology (Steele et al., 2001), and the model is integrated for 1300 years.
There is a rapid adjustment during the first ∼50 years of the simulation, during which the simulated Southern Ocean sea ice extent and upper ocean properties are equilibrating to the surface forcing. After the initial adjustment period, the Weddell Sea water mass properties undergoes a long-term adjustment over the 1300-year simulation period (see Supplementary Figure 1 for the full time series), which is related to the occurrence of two large polynyas forming in the central Weddell Sea. The first polynya event occurs after year 133 in the simulation, persists for 25 years, and is followed by a long period of 300 years without deep convection. The second polynya occurs from year 472 to 500 with approximately the same size and duration as the first event.
In the remaining part of the simulation no polynya forms. Hence, this simulation allows us to investigate how long-term changes in Weddell Sea stratification and WDW heat content influence the formation of Weddell polynyas, and identify key processes leading to the cessation of open-ocean deep convection.

RESULTS
The first polynya initially develops in the deep part of the central Weddell Sea at 23-5 • W and 60-65 • S, hereon referred to as the polynya region (orange box in Figure 1). This region is characterized by a doming of isopycnals, causing a weakening of the local stratification, and possibly promoting the formation of polynyas (Gordon et al., 2007).
In order to diagnose the processes leading to the formation of the polynya in the NorESM-OC1.2 simulation, we limit our analysis to the first 800 years of the simulation. However, the first 50 years are excluded, to avoid any effects associated with the initial model adjustment. First, a detailed analysis of the first polynya event is performed, where we focus on three main time periods: (1) A preconditioning phase (year 50-100) where heat and salt builds up in the WDW layer; (2) A mixing phase (year 100-133) where the WDW starts to erode and sea ice concentration decreases; and (3) A convective phase (year 133-160) marked by the opening of a small polynya which triggers open-ocean deep convection and forcing an abrupt and extensive retreat of sea ice in the Weddell Sea. The sea ice conditions corresponding to each of the different phases are illustrated in Figure 2. After analysing the first polynya event we compare it to the second event, assessing which processes are critical for the polynya to reoccur.
A comparison between the Levitus climatology (Levitus, 1983) and simulated temperature and salinity fields in the Southern Ocean is shown in Supplementary Figure 2. The simulated water mass properties agrees reasonably well with the observational data showing a northward intrusion of low salinity Antarctic Intermediate Water and a southward extension of relatively warm and salty North Atlantic Deep Water (NADW). In the interior Weddell Sea, the WDW is observed as a distinct temperature maximum at mid-depth and originates from entrainment of relatively warm Circumpolar Deep Water (CDW) by the Weddell Gyre (Orsi et al., 1993). However, the simulated WDW layer sits deeper in the water column compared to the observations (see also Robertson et al., 2002;Fahrbach et al., 2004). We attribute this to the relatively weak AABW formation in the model, when deep convection is absent thus resulting in a downward displacement of deep isopycnals causing the core of the WDW to migrate downwards. The implications of the deep WDW layer for the dynamics of the simulated polynya is discussed further in section 4.1. Previous studies have  shown that general circulation models struggle with correctly representing water mass properties in the Southern Ocean and therefore exhibit significant biases relative to observations (e.g., Heuze et al., 2013;Downes et al., 2015). Similarly, in the Arctic, models commonly produce an Atlantic Water layer that is too thick and too deep (e.g., Ilicak et al., 2016). While the simulated stratification in the Weddell Sea region differs somewhat from observations, our model captures the main features of the water mass structure, which allows us to explore the connection between open-ocean polynyas and long-term changes in WDW properties.

Phase I-Preconditioning Phase
The preconditioning phase (year 50-100), corresponds to a period where no polynyas occur. The simulated winter mean sea ice concentration (Figure 2A) and thickness (not shown) in the Weddell Sea is relatively stable throughout the preconditioning phase and compares reasonably well with observations (Parkinson and Comiso, 2008). Thick and compact multi-year sea ice is generally found along the continental shelf and the Antarctic Peninsula due to strong sea ice formation and convergence in the sea ice drift. In the central part of the Weddell Sea, where the polynya initially forms, the seasonal sea ice cover is typically less thick (about 0.6 m) and maximum sea ice concentration ranges between 90 and 100%. The stable winter sea ice cover in the polynya region is sustained by a relatively weak density gradient-the pycnocline-separating warm and salty water in the WDW layer from the surface ( Figure 3C). The upper 200 m, corresponding to the depth of the winter mixed layer, is relative cold and fresh due to freshwater input from seasonal sea ice melt. Below the winter mixed layer, we find the relatively homogeneous and weakly stratified Modified Warm Deep Water (MWDW) which refers to WDW that has been mixed with colder surface waters. The upper limit of WDW, at 1,000 m depth, is characterized by temperatures >0 • C and is separated from the MWDW by a strong gradient in temperature and salinity. Vertical profiles of potential density in Figure 3C illustrate that stratification from the surface down to 2,500 m is strongly dominated by salinity, while decreasing temperatures below the WDW control deep ocean stratification.
During the preconditioning phase we find a gradual buildup of heat and salt in the WDW layer (Figures 4C,D). The accumulation of heat and salt at mid-depth is enabled by a stable stratification, that persists throughout the preconditioning phase and isolates the WDW from surface interactions. At the end of the preconditioning phase (year 100), the core of the WDW layer is located between 1,000 and 2,000 m depth, with a temperature maximum (θ max ) of 0.64 • C and salinity 34.67g kg −1 . This corresponds to a mean warming rate of about 0.005 • C yr −1 and 0.009 g kg −1 yr −1 increase in salinity over the 50 year period. This is within the uncertainties of the observational estimates from Robertson et al. (2002) of 0.012 ± 0.007 • C yr −1 warming of WDW following the 1970s Weddell Polynya. Interestingly, θ max peaks 30 years before the polynya opens, indicating that the build-up of subsurface heat alone is not enough to trigger deep convection. Rather, the heat reservoir preconditions polynya formation by gradually weakening stratification and fuels deep convection once the polynya has formed.

Phase II-Mixing Phase
From year 100, θ max starts to decrease, which marks the onset of the mixing phase. Above the WDW layer, temperature and salinity increases (Figures 3A,B), and the upper ocean stratification weakens (Figure 3C), leading to reduced sea ice concentration in the polynya region. However, the polynya does not appear before year 133 (Supplementary Figure 3).
To understand the processes leading to the opening of the polynya, we look in detail on the time period from year 100 to 150 in Figure 5. The average sea ice thickness in the polynya region is 0.6 m in winter and melts away entirely in the summer (Figure 5A), maintaining a strong surface stratification and facilitating sea ice growth the following winter. Starting in year 125, 8 years prior to the opening of the polynya, sea ice thickness gradually decreases. The reduced sea ice thickness can be linked to an increase in bottom-up fluxes of heat and is suggestive of a vertical diffusion of heat from the warmer FIGURE 4 | Time series of (A) September sea ice concentration (%), (B) total convective area (10 6 km 2 ) and vertical profiles of annual mean (C) potential temperature ( • C) and (D) salinity (g/kg) averaged over the polynya region (23-5 • W, 60-65 • S) for the first 800 years of the simulation. The convective area is defined as the integrated area where the annual maximum mixed layer depth exceeds 2,000 m. In (C,D) the upper 500 m has been exaggerated for a better view of the surface stratification. Black dashed lines mark the opening of the first and second polynya and gray shaded areas show periods of active deep convection.
subsurface layer into the mixed layer. This heat flux melts the sea ice, while maintaining the winter SST at the freezing point ( Figure 5B). When the winter sea ice thickness drops below 0.5 m, the ocean-atmosphere heat flux increases (Figure 5D), thus stimulating further thinning of the sea ice. At the end of the mixing phase, the winter-mean sea ice thickness has decreased to 0.3 m. The annual mean freshwater flux (expressed in terms of sea ice equivalents) associated with melting and freezing of sea ice ( Figure 5E) increases in parallel with the sea ice thinning and increasing heat fluxes. Early in the mixing phase there is a net positive surface freshwater flux, equivalent to 0.3 m of sea ice melt over a year, indicating that more sea ice melts than what is formed locally in the polynya region. Hence, sea ice, that is formed nonlocally, is advected into the region where it melts and freshens the surface layer. As the opening of the polynya approaches, the freshwater flux increases due to a combination of reduced sea ice formation in winter (less brine rejection) and increased advection of sea ice into the region. This tends to stabilize the water column, providing a negative feedback that oppresses further sea ice melt and prevents deep convection (Martin et al., 2013). However, the fact that sea ice thickness continues to decrease shows that the freshwater input is too weak to suppress the strong heat flux from below.
Despite of the increased surface freshwater flux from sea ice melt, there is an increase in the surface salinity prior to the opening of the polynya (Figure 5C). This is explained by the enhanced upward fluxes of warm, saltier water from below, rather than changes in surface fluxes. To demonstrate this, the total salt flux is plotted in Figure 5F and shows contributions from the atmosphere and sea ice model to the surface salinity (i.e., precipitation-evaporation, liquid and frozen run-off, and thermodynamic sea ice growth/decay), but excludes salt fluxes from the ocean model itself. Because the total surface salt flux is negative throughout the mixing phase (mainly due to sea ice melt; see Figure 5E), this indicates that the increase in SSS must be due to an oceanic salt flux into the mixed layer from below and/or horizontally. This preconditions the opening of the polynya in year 133 by weakening the upper ocean stratification (Figure 5G). The mixed layer depth ( Figure 5H) does not change until after the polynya emerges. This supports the fact that the small polynya opening is not triggered by deep convection (through surface buoyancy or wind-stress forcing). Rather, it is only after the polynya has been established that deep convection can occur (Gordon, 1982).

Phase III-Convective Phase
In year 133 a small opening in the sea ice cover appears at 23-5 • W, 60-65 • S during mid-winter (Figure 2; center middle panel). The opening of the small polynya triggers open-ocean deep convection by exposing the underlying ocean to the cold atmosphere, abruptly increasing the ocean-atmosphere heat flux by close to 100 W m −2 . Intense cooling at the air-sea interface leads to enhanced surface buoyancy loss and the water column becomes unstable, triggering deep convection, and mixing over the whole water column ( Figure 3C). As the winter mixed layer deepens (Figure 5H), WDW is injected directly into the surface layer. This leads to winter SSTs well above the surface freezing point (Figure 5B), thus preventing ice formation. A few years after the polynya appears, a large embayment has formed in the Weddell Sea ( Figure 2B; upper right panel).
For an indication of the strength of the simulated deep convection and polynya extent, we calculate total area of convection by summing-up the total number of grid cells where the MLD exceeds 2,000 m following de Lavergne et al. (2014) (Figure 4B). We use the MLD criterion from Martin et al. (2013), corresponding to the depth where the density differs by 0.01 kg m −3 from its surface value. When the polynya opens in year 133, the convective area is only 57,150 km 2 , but grows to a maximum of 2,157,000 km 2 in year 142 ( Figure 4B). In comparison, the 1970's Weddell Polynya was about 300,000 km 2 . This is a common bias in fully-coupled, coarse resolution climate models, producing excessive deep convection and overestimating polynya size (averaging 930,000 km 2 in 70% of pre-industrial CMIP5 model simulations; de Lavergne et al., 2014). We attribute the larger simulated polynya to the substantial ocean heat loss during the convective phase. The change in ocean heat content is shown in Figure 6, and reflects the WDW cooling and diminishing WDW layer thickness. Heat depletion is not only confined to the convection region, but affects the entire Atlantic sector with a heat loss of more than 14 GJ m −2 in the central Weddell Sea. This corresponds to a cooling of WDW by 1.8 • C (θ max drops from 0.6 to −1.2 • C) by the end of the convective phase. The cooling is faster than the resupply of heat from the inflow of CDW, which ceases during the convective phase. In Figure 6A this can be seen as a strong positive anomaly in the South Atlantic corresponding to an accumulation of subsurface heat there. Convection stops in year 160 after the heat reservoir is depleted. Without the upwelling of warm water, the sea ice cover regrows and a cold and fresh surface layer is quickly established thus allowing heat and salt to build up again (Figure 4).
Deep convection associated with the polynya triggers a large increase in bottom water formation (Martinson et al., 1981). This is illustrated in Figure 7, showing the age of water (AOW) at different depths, which is used as a diagnostic for ocean ventilation. A few years before the polynya opens the AOW at 2,000 m decreases, while it increases above the WDW layer at 500 m depth. This points to an increased vertical exchange between the WDW and the more ventilated water mass above. Once the polynya opens we find young water at 4,000 m depth, while near-surface waters are getting older as it is mixed with less ventilated deep water. This a clear indication of open-ocean deep convection which mixes cold, dense surface water down to   Figure 1. When the water is at the surface AOW is reset to zero and older ages are thereby associated with less ventilated water. Vertical dashed black lines shows the opening of the two polynyas in year 133 and 472, respectively. the abyssal ocean. Consequently, the Weddell Sea Bottom Water (below 4,000 m) cool and freshen when convection is active, followed by a mean warming of 0.019 • C per decade, when the polynya is absent (Figures 4C,D). The post-polynya warming in the abyssal Southern Ocean of 0.002-0.019 • C per decade found by Zanowski et al. (2015) is of similar magnitude and is comparable to our result. In response to the stronger deep water formation, the overturning circulation in the Atlantic shows an enhanced northward volume transport of AABW across 30 • S at 3,000 m depth by 7-9 Sv ∼30 years after the polynya occurs (see Supplementary Figure 4). This in good agreement with Zanowski et al. (2015) and Zanowski and Hallberg (2017), who show that it takes 18-25 years for AABW volume transport anomalies, associated with Weddell polynyas, to reach the South Atlantic. Martin et al. (2015) further showed that the enhanced AABW transport associated with Southern Ocean deep convection events may weaken the Atlantic Meridional Overturning Circulation (AMOC) on a multi-centennial time scale. Meanwhile, we only find a small AMOC weakening (2 Sv), indicating that the polynya-driven transport changes are mainly confined to the South Atlantic. This also highlights that openocean polynyas could be an important mode of Southern Ocean ventilation, and implies that the observed warming and shrinking of AABW in recent decades could potentially be linked to an absence of Southern Ocean deep convection (Purkey et al., 2012;Zanowski et al., 2015).

Details of the Trigger Mechanism
In section 3.2, we showed that the opening of the polynya is preceded by increased heat and salt fluxes into the mixed layer leading to warmer and saltier surface waters, weaker stratification and gradual thinning of the winter sea ice cover (Figures 3, 5). This points to diffusive mixing processes as a potential driver for the simulated polynya formation. To illustrate this we show the diffusivity of diapycnal mixing in Figure 8A. The highest values of diapycnal diffusivity (κ = 10 −2 m 2 s −1 ) is generally found in the surface boundary layer, where diapycnal diffusivity is governed by the surface forcing conditions, while in the deep ocean mixing is dominated by the weak background diffusivity (κ = 10 −5 m 2 s −1 ). The strong gradient in mixing between the boundary layer and the interior ocean reflects how stratification acts to prevent warmer subsurface water from entering the surface layer. From year 110 the diapycnal mixing starts to increase at mid-depth, and intensifies with a couple of strong mixing episodes around year 125 and 130. The change in mixing is largest above the WDW and has increased by an order of magnitude from 10 −5 to 10 −4 m 2 s −1 at the end of the mixing phase ( Figure 8B). In the boundary layer, mixing starts to increase 15 years after it increases in the subsurface. This suggest that the enhanced mixing is not driven by changes in surface forcing conditions, but rather arises from destabilizing processes in the deep ocean.
The enhanced vertical mixing results in a gradual increase in temperature and salinity above the WDW interface (Figures 8C,D), while the heat and salt content of the WDW layer is decreasing (Figures 4C,D). This is consistent with the AOW above the WDW layer getting older (Figure 7), while water masses at 2,000 m depth are getting younger. The warmer and saltier water is gradually mixed into the surface layer and temperature and salinity start to increase from year 125. When the polynya opens in year 133, diapycnal mixing is increased over the entire water column mimicking open-ocean deep convection.
As the contribution from shear instabilities to the total diffusivity is generally small (e.g., Dufour et al., 2017) and background diffusivity and tidal mixing is constant, the changes in diapycnal mixing at depth is thought to mainly reflect gravitational instabilities arising from vertical variations in stratification. This is due to the inverse relationship between stratification and vertical mixing (i.e., Equation 2): If the water column is strongly stratified, diapycnal mixing is suppressed, while a weak stratification promotes higher vertical mixing rates. Comparing Figure 8 and Figure 3 illustrates how the parameterization of diapycnal mixing is linked to stratification. The erosion of the WDW layer early in the mixing phase causes a buoyancy flux across the upper WDW boundary that weakens the density gradient and reduces the pycnocline strength (Figures 3C, 5G). This can be attributed mainly to the upward salt flux from the WDW layer, that dominates the stratification at these depths. The background stratification weakens in response, which again favors gravitational instabilities and leads to an erosion of the subsurface heat reservoir. This mechanism presents a positive feedback loop between diapycnal mixing and stratification; a weak stratification favors gravitational instabilities that enhance diapycnal mixing, which in turn weakens stratification further (Dufour et al., 2017).

Description of the Second Polynya Event
The second polynya forms further west in the Weddell Sea near 66 • S, 45 • W (Figure 2B; lower panel). However, the polynya formation mechanism is similar to the first event, e.g., preconditioning-mixing-convection, revealing a multicentennial mode of polynya formation and decay (see also Goosse and Fichefet, 2001;Pedro et al., 2016).
After the first polynya closes in year 160, heat starts accumulating below the surface and the WDW reaches a maximum of 0.1 • C after ∼100 years, similar to the build-up time preceding the first polynya. This corresponds to a mean warming rate of 0.013 • C yr −1 , larger than the warming rate for the first polynya event (0.005 • C yr −1 ), but in agreement with observations (Robertson et al., 2002). This is due to the rapid resumption of CDW inflow resulting in an overshoot in θ max immediately after the first polynya closes (Figure 4C). Similar to the first event, the heat reservoir is maintained by thick and compact sea ice and a pronounced pycnocline which limits vertical exchanges between the WDW and the surface. A gradual deepening of the WDW core, through the preconditioning phase, also contributes to the build-up. We note that the pycnocline strength, indicated by the annual mean salinity difference between the mixed layer and the WDW (Figure 9B), weakens from ∼0.3 psu in year 100 to ∼0.2 psu in year 250. This makes the WDW more vulnerable to surface interactions and could potentially erode the heat reservoir and explain the lower θ max .
Despite the lower WDW heat content it remains warm enough to melt sea ice and preventing ice formation in winter. Thus, when the second polynya appears at year 472 FIGURE 9 | Time series of annual mean (A) potential temperature ( • C) and (B) salinity (g/kg) in the Weddell Sea Convective Region averaged between 31 • W − 9 • E and 60-70 • S at the bottom of the winter mixed layer (black line with squares), the WDW core, e.g., θ max (orange line with triangles) and below the WDW layer at 2,500 m depth (green line with circles). Vertical dashed lines mark the start of each polynya event.
( Figure 2B; lower middle panel) it triggers deep convection and the upwelling of WDW forces a rapid sea ice retreat ( Figure 2B; lower right panel). A few years after the polynya opens, convection reaches its second peak and young surface water is visible at 4,000 m depth, while AOW near the surface gets older (Figure 7). Upwelling ceases in year 496 when the heat reservoir is depleted and sea ice recovers. The heat loss associated with the second polynya event is notably smaller than the first, due to the lower θ max (Figure 6B). Interestingly, the polynya is similar in size and duration as the first polynya and therefore also triggers an increase in AABW of the same magnitude (Supplementary Figure 4). When the second polynya disappears, subsurface heat starts building up again, but only reaches a θ max of −0.6 • C (year 650). Consequently, no polynya forms in the remaining part of the simulation.

Cessation of Deep Convection
To illustrate how the simulated long-term changes in Weddell Sea water mass properties (WDW in particular) and stratification impact polynya occurrence we compare Figure 9 with the vertical buoyancy profiles for each of the preconditioning phases shown in Figure 10. This may help in identifying key processes in polynya formation and understand why deep convection is absent in the later part of the simulation. Figure 10 demonstrates that the stability of the water column (expressed in terms of buoyancy) is controlled mainly by salinity. Here, a fresh surface layer prevents cold water from sinking and the increase in salinity with depth is hindering warm and buoyant subsurface water in rising to the surface. The progressive cooling of the WDW (Figure 9A), therefore acts to stabilize the water column by reducing buoyancy below 1,000 m and weakens the convective potential. Meanwhile, the freshening of the WDW acts to increase buoyancy and weakens the background (haline) stratification at depth, as the salinity difference between the deep ocean (2,500 m) and the WDW diminishes ( Figure 9B). Overall, the balance is in favor of increasing the total buoyancy of the WDW, making the water column less stable. As a consequence gravitational instabilities can be triggered more easily despite the lower θ max value. Above the WDW layer, we find an increase in buoyancy associated with a freshening of the subsurface waters after the first polynya. This tends to stabilize the water column, working against the weaker stratification at depth (Figure 10), FIGURE 10 | Annual mean vertical buoyancy profiles (mm s −2 ) averaged in the Weddell Sea Convective Region (outlined by green box in Figure 1) showing the total buoyancy (black), buoyancy due to the salinity term (blue), and buoyancy due to the temperature term (red) for years 100 (solid) and 250 (dashed). Contributions from temperature and salinity on the total buoyancy are calculated using a linearized equation of state. Here, the buoyancy is expressed as an anomaly with respect to the mean buoyancy profile at the end of each preconditioning phase. The dotted black line shows the total buoyancy at year 650 following the second and final polynya event. and suppress convection for more than 200 years before the second polynya forms.
After the 2nd polynya event, the WDW cools and freshens even further (Figure 9). The water column below 1,000 m is weakly stratified with a salinity difference across the pycnocline of ∼0.15 psu, thus contributing to less heat and salt buildup at mid-depth. Meanwhile, the surface layer (upper 100 m) becomes progressively fresher, which acts to suppress deep convection (de Lavergne et al., 2014). Interestingly, a closer look at Figure 9 shows that temperature and salinity increases above the WDW around year 700, indicative of enhanced vertical mixing. However, at this point the WDW is likely too fresh and too cold to erode the upper ocean stratification and trigger another polynya.

The Role of Mixing in Triggering Polynya Events
Ocean mixing plays a central role in the formation of polynyas simulated by the model, and is strongly linked to long-term changes in WDW properties and stratification. Climate model simulations by Kjellsson et al. (2015) and Timmermann and Beckmann (2004) demonstrate that the Weddell Sea stratification is highly sensitive to the parameterization of vertical mixing, thereby having a direct influence on polynya formation. Insufficient vertical mixing due to inadequate representation of buoyancy-and/or wind-induced mixing can lead to an accumulation of salt in the winter mixed layer that erodes the halocline, causing large and unphysical polynyas to form. Subsequently, it was suggested by Heuzé et al. (2015) that increasing vertical mixing in the upper ocean may reduce deep convection and prevent polynyas from forming. In contrast, we find that the polynya forms in response to enhanced vertical fluxes of heat and salt from below the mixed layer, associated with gravitational instabilities in the ocean interior ( Figure 8A). This implies that increased diapycnal mixing is not related to changes in the surface forcing conditions (e.g., brine rejection, freshwater fluxes, or wind-stress), but rather forced by gradual changes in the WDW properties (Figure 9) thereby increasing the potential instability of the water column. A similar processes is also evident in coupled ocean-atmosphere models where the onset of deep convection is preceded by a gradual erosion of the subsurface heat reservoir over years or decades before the polynya opens (e.g., Martin et al., 2013;Zanowski et al., 2015;Vettoretti and Peltier, 2016;Dufour et al., 2017;Reintges et al., 2017;Zhang et al., 2019). What exactly triggers the convective instability is less clear, but Vettoretti and Peltier (2016) highlighted the role of double-diffusive mixing (not included in the NorESM-OC1.2) as a potential mechanism.
Whether this process plays a role in triggering the observed Weddell Polynya is more uncertain. Based on in situ observations for the 2017 Maud Rise polynya, Mojica et al. (2019) found that diapycnal and isopycnal mixing was important for initiating the polynya, with changes in the subsurface waters occurring before the polynya opened. Here, the enhanced bottom-up fluxes of heat and salt contributed to a higher level of instability in the subsurface waters, thus pointing to an oceanic preconditioning and supports that mixing processes plays an active role in polynya formation. This further highlights the importance of vertical mixing parameterizations to accurately simulate openocean deep convection in climate models (Timmermann and Beckmann, 2004;Heuzé et al., 2015;Kjellsson et al., 2015).
The vertical mixing also plays an important role for the build-up of the subsurface heat reservoir, which in turn affects polynya formation. Throughout the simulation we find a gradual decrease in WDW temperature and salinity (Figure 9) and a cessation of deep convection in the model. This could be linked to weaker stratification and enhanced levels of vertical mixing following each polynya event that prevents heat and salt in building up. This is supported by model simulations by Dufour et al. (2017), who found that heat accumulation relies on a stable stratification where vertical mixing is suppressed. They identified a particular role of representing dense water overflows and mesoscale eddies to simulate the stably stratified water column in the Weddell Sea. When stratification is too weak, arising from a poor representation of water masses, episodes of gravitational instability can trigger enhanced vertical mixing and lead to erosion of the WDW reservoir.
Another factor that may affect the subsurface heat reservoir is the deepening of the WDW layer. A shallower WDW core, is more vulnerable to upper ocean mixing, leading to increased iceocean heat fluxes that deplete the heat reservoir (e.g., Robertson et al., 2002). As evident in Figure 4, the simulated WDW is located around 1,000 m depth, and is therefore less influenced by surface forcing. We attribute the WDW deepening to a weak resupply of bottom water during the non-convective phases. When the polynya and deep convection is absent bottom water formation is reduced and we find a downward migration of the WDW core at a rate of 12 m yr −1 (Figure 4). This helps isolate the WDW and may also create favorable conditions for the formation of a polynya through cabelling and thermobaric effects (Akitomo, 1999;Mcphee, 2003;Wang et al., 2016). A similar response has been documented in other higher resolution models (e.g., Lee et al., 2002;Dufour et al., 2017) and represents an ongoing challenge influencing longterm climate model predictions. For example, Dufour et al. (2017) found an isopycnal displacement of 10 m yr −1 which was attributed to a lack of resupply of AABW and erosion of AABW properties. Alternatively, isopycnals models have been shown to be sensitive to the choice of the reference pressure, in particular in the Southern Ocean, (e.g., Assmann et al., 2010). Hence, the choice of 1,000 m as a reference level for our model density layers may cause stratification biases in the abyssal ocean and enhance consumption of bottom waters by spurious mixing processes.

What Controls the Duration and Frequency of Polynya Events?
We find two time scales associated with the frequency of openocean deep convection in NorESM-OC1.2. First, an advectivediffusive time scale associated with accumulation of subsurface heat, and second, a time scale set by the stratification and ocean mixing processes. The build-up of the heat reservoir occurs in about 100 years and is associated with advection of CDW into the Weddell Sea. Martin et al. (2013) show a similar time scale in the Kiel Climate Model and found that the centennial-scale recharge of the subsurface heat reservoir can drive an oscillatory mode of Southern Ocean deep convection, with the heat reservoir setting the time between polynya events. Meanwhile, observations suggest that already by the mid-1990s the WDW had recovered from the cooling induced by the 1970s convective event, which is significantly faster than the 100 years in our simulations (Robertson et al., 2002;Smedsrud, 2005). While observations remain too sparse to asses the full spatial extent of the WDW cooling by the Weddell Polynya, hydrographic profiles suggest that a significant amount of heat is left in the WDW. This lead to a faster recovery of the heat reservoir and the recurrence time for a Weddell Sea polynya is expected to be faster, than if the entire heat reservoir is depleted. Due to the large size of the simulated polynya, the total oceanic heat loss associated with convection is about 3.6 × 10 22 J. This is an order of magnitude larger than the heat loss from 1974 to 1976 during the Weddell Polynya (Smedsrud, 2005) and is likely due to a combination of poor representation of subgridscale convection in coarse resolution models and the lack of atmospheric feedbacks. Despite the large simulated oceanatmosphere heat fluxes associated with deep convection, the atmosphere remains cold and the strong surface buoyancy loss continues to drive upwelling of WDW until the heat reservoir is depleted. Including atmospheric feedbacks could help dampen convection by restratifying the water column, e.g., through changes in wind-driven upwelling, heat and moisture fluxes, or sea-ice transport into the polynya region (Martin et al., 2013;Gnanadesikan et al., 2020). Therefore, the 100-year time scale for the recharge process should be seen as an upper limit. Nevertheless, this result suggest a close link between the WDW heat content, polynya size and the frequency of polynya events: Large and long-lived polynyas leave the heat reservoir more depleted, the recovery time longer, and deep convection occurs less frequent. In contrast, for short-lived and smaller polynyas less heat is lost and the recharge time is shorter.
Changes in the local atmospheric conditions and oceanatmosphere feedbacks are other process that could influence deep convection activity by affecting Weddell Sea stratification (e.g., Gordon et al., 2007;Gnanadesikan et al., 2020;Kaufman et al., 2020). For example, surface freshening by wind-driven advection of sea ice, or precipitation anomalies, acts to stabilize the water column and suppress polynya formation (Comiso and Gordon, 1987;de Lavergne et al., 2014). On the other hand, recent analysis of the 2016 and 2017 polynya event, indicate that variations in Southern Hemisphere winds were important for triggering the polynya, by creating ice divergence and wind-driven upwelling of WDW (e.g., Cheon et al., 2015;Campbell et al., 2019;Francis et al., 2019). Such variations may be linked to long-term changes in large-scale atmospheric patterns like the Southern Annular Mode (Gordon et al., 2007;Jena et al., 2019). However, since the atmospheric forcing is fixed in our simulation, there is no external variability in the surface forcing conditions (i.e., winds, precipitation) that could impact polynya activity. Instead, changes in stratification are driven by internal ocean processes (see also Dufour et al., 2017) and thus presents a purely oceanic time scale for polynya formation. This absence of stochastic forcing by the atmosphere likely leads to polynyas forming less frequently. Therefore, we speculate that the oceanic processes described here are dominating polynya activity on centennial time scales (through advective-diffusive processes) whereas atmospheric variability is important over shorter inter-annual to decadal time scales.

Concluding Remarks
Using a stand-alone ocean-sea ice model (NorESM-OC1.2), we show that internal ocean processes can trigger the formation of open-ocean polynyas and drive deep convection in the Weddell Sea, without changes in the atmospheric forcing. The proposed mechanism relies on the build-up of WDW as an essential preconditioning for polynya formation (Cheon et al., 2015;Dufour et al., 2017), leading to gravitational instabilities in the ocean interior and enhanced vertical mixing.
The polynya formation mechanism can be separated into three phases: (1) a preconditioning phase characterized by accumulation of WDW at mid-depths sustained by a strong stratification that suppresses vertical mixing; (2) a mixing phase associated with enhanced upward fluxes of heat and salt leading to a weakening of the upper-ocean stability and a small polynya opening; and (3) a convective phase, where deep convection drives upwelling of WDW that maintains the polynya. After a couple of decades, the subsurface heat reservoir is depleted and the system returns to phase (1).
Finally, we show that deep convection causes significant changes in Weddell Sea deep and bottom water properties which eventually leads to a shut-down of deep convection in the model. This provides a perspective on the long-term evolution of Southern Ocean polynya activity. We find that the polynya cannot form when the WDW becomes too cold (θ max <−0.6 • C). In comparison, Cheon et al. (2014) found that open-ocean polynyas can form even for very low θ max values close to the freezing point (−1.4 • C), implying that additional factors (such as the salinity of WDW) might be important for initiating the polynya. Here we show that the destratification by subsurface ocean mixing is prevented when the WDW is freshened, thereby suppressing ocean-ice heat fluxes and sea ice melt. Therefore, we propose that the combination of a colder and fresher WDW may lead to a cessation of deep convection. As opposed to the changes simulated here, the past decades have shown a warming of the Southern Ocean deep waters (e.g., Purkey and Johnson, 2010;Fahrbach et al., 2011), which could help destabilize the water column and favor deep convection. At the same time global climate models predict convection to decrease in the twenty-first century due to a continued freshening of Southern Ocean surface waters under warming conditions (de Lavergne et al., 2014). Given that the warming trend of WDW is enough to overcome the enhanced surface stratification, the likelihood of large and long-lived Weddell polynyas occurring in the future could therefore increase.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
The analysis, plotting and preparation of the manuscript was done by JR. LS and KN also aided with the interpretation of the model simulation. All authors contributed to the writing and editing of the final manuscript.

FUNDING
The research leading to these results has received funding from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013)/ERC Grant agreement 610055 as part of the ice2ice project.