Causes and Evolution of Winter Polynyas over North of Greenland

During the 42-year period (1979-2020) of satellite measurements, only three winter polynyas 15 have ever been observed north of Greenland and they all occurred in the last decade, i.e. February of 2011, 2017 and 2018. The 2018 polynya was unparalleled by its magnitude and duration compared to the two previous events. Combined with the limited weather station and remotely-sensed sea ice data, a fully-coupled Regional Arctic System Model (RASM) hindcast simulation was utilized to examine the causality and evolution of these recent extreme events. We found that neither the accompanying 20 anomalous warm surface air intrusion nor the ocean below had an impact on the development of these winter open water episodes in the study region (i.e., no significant ice melting). Instead, the extreme atmospheric wind forcing resulted in greater sea ice deformation and transport offshore, accounting for the majority of sea ice loss. Our analysis suggests that strong southerly winds (i.e., northward wind with speeds of greater than 10 m/s) blowing persistently for at least 2 days or more, were required over the 25 study region to mechanically redistribute some of the thickest sea ice out of the region and thus to create open water areas (a latent heat polynya). In order to assess the role of internal variability versus external forcing of such events, we additionally simulated and examined results from two RASM ensembles forced with output from the Community Earth System Model (CESM) Decadal Prediction Large Ensemble (DPLE) simulations. Out of 100 winters in each of the two ensembles, initialized 30 years 30 apart, one in December 1985 and another in December 2015, respectively, 17 and 14 winter polynyas were produced over north of Greenland. The frequency of polynya occurrence and no apparent sensitivity to the initial sea ice thickness in the study area point to internal variability of atmospheric forcing as a dominant cause of winter polynyas north of Greenland. We assert that dynamical downscaling using a high-resolution regional climate model offers a robust tool for process-level 35 examination in space and time, synthesis with limited observations and probabilistic forecast of Arctic events, such as the ones being investigated here and elsewhere. https://doi.org/10.5194/tc-2021-279 Preprint. Discussion started: 29 September 2021 c © Author(s) 2021. CC BY 4.0 License.


Introduction
The Arctic has experienced amplified warming, both through the enhancement of global temperature rise, as well as through the reductions in sea ice and snow cover that impact the regional energy budget 40 (Serreze and Francis, 2006). This is commonly referred to as Arctic amplification (AA). On a seasonal basis, Arctic winter warming (AWW) exceeds summer warming by about a factor of four (Bintanja and van der Linden, 2013). In addition to sea ice variations, changes in the atmospheric circulation have been linked to an increase in frequency and duration of winter warming events in the Arctic . The trends in AWW and winter sea ice extent (SIE) have continued over the satellite record, 45 with the March SIE decline rate of ~2.6% per decade for the period of 1979-2020, including the four lowest winter maxima of SIE in 2015-2018 and the lowest ten during 2005-2019. While SIE reductions in winter are less than ones in summer (Stroeve and Notz, 2018), there is clear supporting evidence that the winter ice pack has thinned and the amount of multi-year sea ice has declined (Kwok et al., 2009;Meier et al., 2014;Kwok, 2018). It is expected that the resulting younger and thinner ice is more 50 susceptible to atmospheric wind forcing (Spreen et al., 2011;Itkin et al., 2017), yielding increased sea ice drift speed, enhanced fracturing and more lead openings (Rampal et al., 2009). Hence, along with the current trend in the Arctic sea ice toward younger and thinner ice, it can be hypothesized that polynyas may become more prevalent in recent years. Since about half of the total atmosphere-ocean heat exchange over the Arctic Ocean may occur through leads and polynyas in winter (Maykut, 1982), 55 the occurrence and formation of polynyas could play a crucial role in the alteration of regional climate (Morales Maqueda et al., 2004).
In late February 2018, satellite imagery revealed an unusual open water polynya over north of Greenland between the Lincoln and Wandel seas (Fig. 1a). This event received considerable attention 60 not only because it was claimed to be a one-of-a-kind extreme event involving some of the thickest Arctic sea ice but also because its emergence coincided with anomalous, above freezing, warming of surface air temperature over the region (Moore et al., 2018). The presumed contribution of this warm surface air to the polynya opening is of interest and so are its causality, evolution and past occurrences. Yet, detailed in situ observations of polynyas are limited due to their intermittency and restricted access 65 in winter. Hence, in addition to satellite measurements and weather data, fully coupled climate models become critical tools in studying such events (e.g., Ludwig et al., 2019). However, for the majority of global climate models (GCMs), many coastal polynyas are at sub-grid scales; thus, the GCM utility for comprehensive polynya studies are impeded (Weijer et al., 2017). In addition, GCMs are not intended to represent specific climate events in space and time and they are not suitable for process-level 70 investigation of such events and quantification of their impact at both local and larger scales.
On the other hand, regional climate models (RCMs) used for dynamical downscaling are expected to show improvement in reproducing extreme weather events, given that their atmospheric boundary conditions for simulations of the past to present are derived from global atmospheric reanalysis such as 75 Climate Forecast System (CFS) version 2 (CFSv2; Saha et al., 2011). Therefore, high-resolution RCMs, such as the Regional Arctic System Model (RASM; e.g. Maslowski et al. 2012), offer unique capabilities for examining the spatio-temporal development and impact of observed specific events like a polynya (Fig. 1b), in the context of a fully-coupled climate system model (atmosphere-sea ice-oceanland), while CFSv2 is shown to be less skillful in simulating such an event (Fig. 1c). In addition, RCMs 80 afford ensemble sizes prohibitive to their global fine-resolution counterparts, which is often a requirement to distinguish the forced response from internal model variability (e.g., Peings et al. 2021).
The 2018 winter polynya over north of Greenland has been investigated by Moore et al. (2018) using an ice-ocean model forced by surface atmospheric reanalysis. They have established the dominant role of 85 surface winds in generating this polynya. Here, by taking advantage of the fully-coupled and highresolution RASM hindcast simulation, combined with weather station and satellite sea ice data, we evaluate the capability of RCM in reproducing the observed natural phenomena (Fig. 2) and investigate the coupled mechanisms involved in the development of northern Greenland coastal polynyas within some of the thickest Arctic ice-pack cover. In particular, this study focuses on analysis of historical 90 winter polynya events for the full period of satellite data availability  to diagnose the relative roles of thermodynamic and dynamic processes and to assess required forcing changes over the last four decades. Furthermore, by dynamic downscaling of the Community Earth System Model (CESM)-Decadal Prediction Large Ensemble (DPLE), two sets of 30 years apart RASM ensemble (termed RASM-DPLE) simulations are performed to examine the relative roles of internal variability 95 and external forcing influencing the development of winter polynyas under two different sea ice regime scenarios: i.e. a thicker ice in the 1980s versus a thinner ice in the 2010s. We provide details of the satellite and weather station data used for this study in section 2 and the model setup for the hindcast and DPLE simulations are described in section 3. Next, section 4 presents a synthesis of observed and modeled past winter polynya results and examines the statistics of polynya occurrence and the required 100 conditions for their generation in the RASM-DPLE simulations. This is followed by the discussion in section 5 and the study is summarized in section 6. Daily sea ice concentration (SIC) data were obtained via the National Snow and Ice Data Center (NSIDC): https://nsidc.org/data/nsidc-0051 for 1979-2019 and https://nsidc.org/data/nsidc-0081for 2020. Satellite-derived SIC used for this study is based on passive microwave measurements using the NASA team (NT) algorithm (Cavalieri et al., 1984); all the data are on a polar stereographic 25 km×25 120 km grid. The mean daily SIC was used to detect the occurrence of a polynya in the satellite measurement when it dropped below 90% over the study region. RASM sea ice thickness (SIT) was also compared with the CryoSat-2/the Soil Moisture Ocean Salinity (SMOS) satellite merged data that are only available in winter months from 2011 (https://www.meereisportal.de; Grosfeld et al., 2016) as well as CFSv2 reanalysis (https://doi.org/10.5065/D61C1TXF). Due to the lack of persistent SIT 125 observations over the Arctic, the Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS) is often considered as an "observational" proxy (Zhang and Rothrock 2003;Schweiger et al. 2011;Stroeve et al. 2014 Maslowski et al. 2012;Roberts et al. 2015;DuVivier et al. 2015;Hamman et al. 2016;Hamman et al. 2017;Cassano et al. 2017). All the components are coupled using the Craig et al. (2012) version of CESM flux coupler. RASM is configured over a pan-Arctic domain, including the entire 140 Northern Hemisphere marine cryosphere and all terrestrial drainage basins that drain to the Arctic Ocean. The ocean and sea ice components share a 1/12-degree (~9 km) rotated sphere grid and are configured with 45 levels in the vertical and five thickness ice categories, respectively. The atmosphere and land hydrology components are set up on a 50-km polar stereographic mesh with the vertical resolution of 40 levels and 3 soil layers, respectively. The RASM-DPLE simulations are derived by dynamically downscaling global atmospheric output from 155 the initialized CESM-DPLE simulations (Yeager et al., 2018). Output from the two 10-member decadal ensembles, initialized on December 1 st of 1985 and 2015, was selected for in-depth analysis under different regimes of the Arctic SIT distribution. Each RASM-DPLE simulation is initialized with the ocean and sea ice conditions, with thinner ice in the latter period (Fig. S1), from the RASM hindcast and integrated for 121 months with CESM-DPLE atmospheric forcing. The size of each ensemble (10 160 members) is determined by the availability of archived CESM-DPLE output necessary for RASM-WRF boundary conditions. Hence, output for 100 winters per each ensemble allows for statistical analysis of polynya occurrence in the past and near future.

Self-organizing maps
The self-organizing map (SOM), an artificial neural network based on a competitive learning algorithm 165 (Kohonen, 2001), has been widely used to visualize input data vectors onto a low dimensional map of nodes and to objectively classify complex data sets in meteorology and oceanography (see Liu and Weisberg, 2011). For January-March 2018, the RASM/WRF pan-Arctic (>65 °N) 6-hourly surface wind fields (i.e., 10 m U and V components) are characterized using a SOM [4×4] map grid (i.e., 16 nodes or patterns). Time-dependent spatial features of near surface winds are identified and frequencies 170 of occurrence of wind patterns favorable for a polynya are quantified using the SOM Toolbox 2.0 for MATLAB available at http://www.cis.hut.fi/projects/somtoolbox.

Near-surface air temperature around Greenland 175
Focusing on the most recent and the largest winter polynya event, daily near-surface air temperatures were examined from the weather stations around the Greenland coast for January-March 2018 (Fig. 3). A significant warming event (air temperatures rising above 0 o C) was observed over the northeastern Greenland region from mid-February to early March and captured well in the RASM simulation. This anomalous warming was most prominent along the northern Greenland coast (Figs. 3e and 3f), where 180 the polynya was observed and simulated in RASM (Figs. 2a and 2d, respectively), but less pronounced over the mid-eastern Greenland coast (Figs. 3g and 3h). In contrast, no warming was measured in the southwestern stations for the same period (Figs. 3b-3d). This anomalous warming, with relative humidity rising above 90% (i.e., Station 04312; not shown), coincided with a strong reversal of the Arctic Oscillation (AO) index from a positive to a strongly negative phase (Fig. 3a). Among the 185 northeastern stations, the observed near-surface air temperature was positively correlated to the AO index (Fig. 3i). The maximum correlation coefficient (r) was time-lagged up to 11 days at the northern stations and only 3-4 days at the mid-eastern stations (Fig. 3j). Given that the anomalous warming started a few days earlier at the mid-eastern stations and ceased a few days later at the northern stations, it may suggest the advection of warm air masses from the south. The prolonged warming at the northern 190 stations could also be partly due to the release of oceanic heat from the polynya to the atmosphere.
On the other hand, near-surface air temperatures were inversely correlated with the AO index ( Fig. 3i) in the southwestern Greenland region, with a shorter time-lag (one to two days) for a maximum correlation coefficient (Fig. 3j). We found that no anomalous warming was present in the southwestern 195 Greenland region during February. In fact, near-surface air temperatures tended to gradually decrease from January through February, with the lowest temperatures observed between 22 nd and 24 th of February 2018 (Figs. 3b-3d). Thereafter, near-surface air temperatures rapidly increased at all southwestern stations, peaking in early March, corresponding to the strong negative AO phase, and the second warming came approximately two to three weeks later. 200 Overall, the RASM hindcast simulation captured remarkably well the sudden increase of near-surface air temperatures in northeastern Greenland and its gradual decrease. RASM also reproduced well the cooling in February and then the warming over southwestern Greenland, in terms of its magnitude, timing, and spatio-temporal variability (Fig. 3). One discrepancy in the RASM simulation against the 205 weather stations data was a positive bias of near-surface air temperatures at Stations 04221 ( Fig. 3b) and 04330 ( Fig. 3g), possibly linked to the relatively coarse horizontal resolution of the RASM atmospheric component (i.e., 50 km), which is insufficient for resolving strong temperature gradients across the ocean/land/ice sheet boundary or the fidelity of the near surface temperature distribution over the Greenland Ice Sheet. 210

Sea ice dynamics
The sudden anomalous warming over northern Greenland with temperatures above 0 o C was an extreme phenomenon in 2018, considering that the long-term mean (2011-2020) of February surface air temperature is -27.8 o C at Station 04301 ( Fig. 3e; Cape Morris Jesup). In addition, there were other years of anomalous warming (Figs. 4b and 4c), albeit less pronounced, when previous, smaller, polynya 215 events occurred during February 2011 and2017 (Figs. 2b and2c, respectively). The RASM's realistic representation of the polynya, although the center of the simulated polynya is not exactly located as the observed one, as well as the magnitude and timing of anomalous warming, grants confidence to the examination of relative contribution of thermodynamic ice melt to the generation of the polynya. Figure  5 shows the RASM thermal sea ice surface, lateral and bottom melting terms were all negligible (< 1 220 cm) over the study region when integrated for the whole month of February 2018. Hence, the dramatic rise of near-surface air temperatures by more than 25 o C above climatology and their persistence around the freezing point for several days (Fig. 3e), had no impact on sea ice melt, nor on the preconditioning and development of the polynya north of Greenland. In agreement with Moore et al. (2018), we corroborate that this polynya was driven by mechanical redistribution of sea ice outside of the study 225 region (see Fig. 4a): i.e., this was a latent heat polynya. Based on RASM results, we calculate that between 15 and 25 February 2018, 192 km 3 of sea ice was dynamically transported outside the study region (Fig. 4d). During two weeks prior to the 2018 polynya event, a mean thermodynamic ice growth over the region was 0.72 km 3 /day (Fig. 4d), which is comparable to the rate in a non-polynya year such as 2019 (not shown). The peak sea ice growth of 3.1 km 3 /day occurred on 26 February right after the 230 maximum daily dynamic ice removal and the anomalous warming period. However, this large sea ice removal due to the polynya was not fully replenished in the region by the end of March, even with dynamic and thermodynamic processes adding 81 km 3 and 42 km 3 , respectively, after 26 February.
Overall, the RASM integrated thermodynamic ice growth in the study region during February 2018 was approximately 50% higher (31.2 km 3 ; Fig. 4d) due to the rapid ice growth during the polynya opening, 235 compared to the ice growth during a non-polynya year (i.e., 20.8 km 3 in February 2019; not shown).

Atmospheric-sea ice coupling
As dynamic processes dominated the overall winter sea ice in the study region, we have analyzed the wind data from the weather station (Station 04301; Station Nord) and have found that the polynya development was associated with strong and persistent winds from the south-southeast (Fig. 4g). We 240 further examined how the spatial near-surface wind fields evolved over the time period associated with the sea ice divergence. The SOM analysis extracted 16 patterns from the total of 360 synoptic wind fields: the 6-hourly RASM 10 m U-and V-wind components during January-March 2018. Figure 6 shows the four major wind patterns most frequently identified (more than 77% occurrence) from the pre-polynya period (5 February  RASM hindcast simulation confirmed that this polynya event was predominantly associated with southerly to southeasterly winds blowing over the northern Greenland region ( Fig. 6b and 6c), consistent with the ERA-Interim reanalysis of 10 m wind fields (Fig. S2). Prior to the polynya event in 2018, surface winds were mainly from the north or northwest over the region (Fig. 6a), which yielded little or no sea ice divergence (Fig. 6e). When the major wind pattern shifted to between southerly to 250 southeasterly winds over the northern Greenland region on 15 February 2018 (Fig. 6b), sea ice started to deform and diverge significantly (Fig. 6f). Beginning on 20 February 2018, the southeasterly wind became even more prominent, with 19% stronger wind speed over the northern coast of Greenland (Fig.  6c), which increased the sea ice deformation rate further and led to the maximum polynya opening (Fig.  6g). The largest observed and modeled polynya areas were identified on 25 February 2018 (see Figs. 2a  255 and 2d, respectively). Thereafter, within a week, the shift of wind patterns in late February (Fig. 6d) reversed the dynamic sea ice volume (SIV) tendency from net loss to gain (Fig. 4d) and subsequently reduced the deformation rate back to nearly zero by early March (Fig. 6h)

Winter polynya events before 2018 260
Upon examining satellite-derived SIC over the entire data record between 1979 and 2020 (Figs. 7, S3, and S4), we found that there were two additional polynya events, albeit smaller, over the northern Greenland region observed in February 2017 and2011 (see Figs. 2b and2c, respectively). Here, we found that a polynya event occurred when the daily mean satellite-derived SIC fell down to or below 90%, which occurred for three consecutive days in 2011, one day in 2017, and 10 days in 2018 (see 265 Figs. 7a, 7g, and 7h, respectively). As was the case of the 2018 polynya, both the previous events also coincided with anomalous warming, peaking on 12 February 2011 (Fig. 4b) and on 8 February 2017 (Fig. 4c), as measured at Station 04301 (Cape Morris Jesup). Near-surface air temperature variability was statistically correlated with the AO index (r=0. 39, p<0.01 in 2011; r=0.45, p<0.01 in 2017), lagged by approximately two weeks. Those two polynyas were represented well in the RASM hindcast 270 simulation (Figs. 2e and 2f) and thus we investigated their causality with respect to their relative SIV reductions due to dynamical processes and/or thermodynamic ice melt (Figs. 4e and 4f). The RASM results confirmed that those were latent heat polynyas dominated by the mechanical redistribution of sea ice out of the region. Their sizes were smaller compared to the 2018 polynya, with the one in 2017 even smaller than the one in 2011, which we attribute to somewhat different wind patterns such as its 275 direction, magnitude and duration ( Figs. 4i and 4h). Table 1 shows that, during the polynya periods (defined here when ice volume tendency is negative), i.e., 12-15 February 2011 and 8-10 February 2017, 55 km 3 and 42 km 3 of sea ice was dynamically removed outside the study region, respectively, which is much less compared to the ice loss of 189 km 3 during the 2018 event (16-25 February 2018). The size of a polynya was proportional to the pseudo-wind stress (i.e., wind speed squared) integrated 280 over the polynya period (Table 1). Thus, more turbulent (latent plus sensible) heat was lost during the 2011 winter polynya (daily mean of -96.7 W/m 2 and maximum of -125 W/m 2 ) when its size was bigger (Table 2). In addition, referenced to the February ice growth during a non-polynya year (for example, 20.8 km 3 in February 2019; not shown), the RASM thermodynamic ice growth integrated over the month of February was elevated: 33% higher in 2011 (27.6 km 3 ; Fig. 4f) in the study region. Note that 285 the daily mean turbulent heat flux in the study region was -48 W/m 2 during the 2018 winter polynya event with the maximum daily heat loss up to -182 W/m 2 . However, the total daily turbulent heat loss was much larger than any other year because of the size and duration of open water areas (Table 2).
Given the above analysis, an outstanding question is why winter polynyas was only observed during the 290 2010s within the past four decades. Since observational data around northern Greenland are incomplete for the whole period of 1980-2020, we expanded the analysis of surface wind fields (10 m U and V components) from the RASM hindcast simulation near Station 04312 ( Fig. 8a; Station Nord) and Station 04301 ( Fig. 8b; Cape Morris Jesup). Our analysis revealed that northward wind (blowing from south), required for opening of a winter polynya along the coast of northern Greenland, has recently 295 become more frequent, stronger (i.e., wind speed >10 m/s), and more persistent (i.e., blowing for at least 2 consecutive days or longer). In addition, the three years of winter polynya occurrence with the wind conditions satisfying the above criteria, all coincidently were La Niña winters (see https://origin.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ONI_v5.php). Based on the RASM hindcast simulation, the grid mean wind conditions in February 2009 were similar to the ones in 300 February 2017 near Station 04312 (Fig. 8a), but a notable polynya was not detected, possibly owing to the influence by such wind conditions over a smaller area within the main polynya region. The observational data also indicate a polynya favorable wind condition in 2009: i.e., relatively warmer air (-12.6 °C) blowing from the south-southwest with the maximum wind speed of 19 m/s on February 7 th , 2009 at Station 04312 (Table S1). But, the satellite-derived mean daily SIC only dropped down to 94% 305 in February 2009 (Fig. S4o) because the wind was possibly weaker, compared to the wind condition in February 2017 (Fig. 4i)  Taking as the baseline the winter polynya in February 2017 ( Fig. 4e and Table 1), which was the smallest among the three observed events, we defined "a latent heat polynya" in the RASM-DPLE ensemble members when a daily winter sea ice loss due to dynamic processes was greater than 10 335 km 3 /day for at least three consecutive days over the study area (see Fig. 4a). Note that the two other observed polynya events experienced 4 and 10 consecutive days of dynamic sea ice loss greater than 10 km 3 /day during February of 2011 (Fig. 4f) and 2018 (Fig. 4d), respectively (Table 1). Analogous to the RASM hindcast simulation, Table 3 and 4 list all the polynya occurrences from the two ensemble runs. We found 17 polynyas in the 1985-initialized ensemble (Table 3) and 14 polynyas in the 2015-initalized 340 ensemble (Table 4), out of the 100 winters each. Note that some ensemble members simulate more than 1 polynya event while some had no polynya at all. The majority of the polynya events (88% in the 1985-initialized and 64% in the 2015-initialized runs) were similar in size to the smaller polynyas that were simulated in February of 2011 and 2017. The longer-lasting polynyas (with dynamic sea ice removal of 10 km 3 /day or greater for five consecutive days or longer) were also produced in the RASM-345 DPLE simulations: two incidents in the 1985-initialized and four incidents in the 2015-initialized runs. But none of the ensemble members produced the total integrated turbulent heat flux as large as the 2018 polynya (Table 2).
We finish this analysis of RASM-DPLE results by comparing the largest polynya detected in each 350 ensemble (see Table 3 and 4), i.e., from the ensemble member #4 in January 1988 (Fig. 9a)  ensemble member #2 in January 2024 (Fig. 9b). As in the case of observed events, these latent heat polynyas were created due to dynamic sea ice transport (DVT ≤ -10 km 3 /day for seven consecutive days), resulting in significant sea ice removal of -120 km 3 out of the region in 1988 ( Fig. 9c; Table 3) and -136 km 3 in 2024 ( Fig. 9d; Table 4). By cross-examining wind patterns over northern Greenland 355 (i.e., near Cape Morris Jesup), our analysis confirmed that these large polynya openings were also associated with very strong and persistent southerly winds in the RASM-DPLE simulations (Figs. 9e and 9f). We found that daily DVT and meridional wind were significantly correlated, and the size of polynyas was highly dependent on the pseudo-wind stress integrated over the polynya period (Table 3 and 4). However, given the wind patterns, their duration and the integrated sea ice removal (-189 km 3 ) 360 during the observed polynya in February 2018 (Table 1), this event stands out within all the RASM results analyzed in this study. Compared to the 2018 event, the slight difference in the magnitudes between the largest polynyas of each ensemble still does not imply much significance of SIT in their generation and evolution. Overall, a similar number of winter polynyas, produced between the two 30year apart ensembles, implies that changes in SIT are not significant contributors (at least up to now) to 365 the generation of such events for this region during wintertime.

Discussion
The dynamic downscaling approach using the high-resolution RASM allows resolving the fine-scale processes and provides valuable insight into the mechanism of generation and evolution of winter polynyas off the northern coast of Greenland. We examined the occurrence of such events in RASM 370 against the satellite observations over the past four decades. The results from the RASM hindcast simulation suggest that the size of a winter latent heat polynya in this region is sensitive to the direction, magnitude and duration of near-surface winds. Subject to the limited sample sizes, we find that the generation of a winter polynya in this region requires strong southerly winds (i.e., speeds greater than at least 10 m/s and lasting for more than 2 days largely over the study region based on the RASM hindcast 375 simulation). Table 1 and Figure 8 show that the stronger and more persistent the southerly wind blows, the larger the winter polynya becomes; however, there were no such extreme weather events before 2010. On the other hand, the southerly winds might reduce sea ice export through Fram Strait by a slowdown of sea ice drift (Wang et al., 2021). It can also be speculated that this region might have been affected by increasing extreme winter storm activities in recent years that are associated with anomalous 380 warming events. However, no significant trends were found in January- February during 1979-2015.
Considering the overall vulnerability and fate of some of the thickest sea ice in the Arctic under the recent warming climate, sea ice may become susceptible to modulation by the atmospheric forcing. 385 However, the RASM-DPLE simulations indicate that the frequency, size and duration of winter polynyas is not affected by sea ice thinning, showing no apparent difference in such polynya characteristics between the two ensembles initialized with two different sea ice regimes 30 years apart: 1980s vs 2010s. For example, the range of 5-day mean SIT before the occurrence of each polynya was 2.9-4.4 m (with the mean of 3.8 m and standard deviation of 0.40 m) in the 1985-initialized RASM-390 DPLE ensemble and 1.8-3.3 m (with the mean of 2.8 m and standard deviation of 0.35 m) in the 2015initialized RASM-DPLE ensemble (Table 3 and 4, respectively). Hence, regardless of the SIT decline over the region, this study suggests that the primary necessary condition for a winter polynya occurrence is strong and persistent southerly winds. With the maximum wind speed exceeding 20 m/s and duration for 10 days, the 2018 winter polynya remains unique by any of these metrics. 395 During the 2017-2018 winter, the AO shifted from a positive phase to a strong negative phase between mid-February to early March (Fig. 3a), which was associated with weakening of the polar vortex and allowed warmer air into the Arctic. Alternatively, Kim et al. (2014) argued that due to sea ice loss especially in Barents-Kara seas, the weakened stratospheric polar vortex preferentially induced a 400 negative phase of the AO at the surface, resulting in warm air into the Arctic. This anomalous warming event occurred coincidently after the sudden stratospheric warming (SSW) event observed in mid-February 2018 (Moore et al., 2018;Rao et al. 2018), which developed into a vortex split (Lü et al., 2020). Subsequently, the winter weather was severe with intense cold air across Europe in March 2018 (Overland et al., 2020). Although the exact cause of SSW variability in still under debate, SSWs are 405 generally known to cause anomalous warming over Greenland and impact surface weather patterns down to mid-latitudes . The frequency of their occurrence is enhanced by El Niño conditions (Polvani et al., 2017 and also true for La Niña winters depending on SSW definitions (Song and Son, 2018); for example, La Niña was in winter 2017-2018. In the recent winters with SSW Rao et al., 2018;Knight et al., 2021), no winter polynya events occurred except in February  410 2018.
The opening of polynyas in the region between the Wandel and Lincoln seas primarily depends on a large-scale surface pressure pattern change, resulting in strengthening of intense southerly winds that are short-lived and sporadic. Toward the end of each polynya event, the rate of thermodynamic ice 415 growth generally increased (Figs. 4d, 4e and 4f). The maximum rate of ice growth happened not when sea ice removal was largest with strong southerly wind but when air temperature started to significantly drop (Fig. 4). Upward surface turbulent heat fluxes continued even after DVT became positive (i.e., Table 2 and Fig. 4) when wind direction was reversed (i.e., Fig. 6d) until the open water area was completely covered by ice. In general, the rate of new ice formation in winter is expected to be high 420 during an open water phase, such as leads and polynyas, due to an intense turbulent heat (sensible and latent) loss to the atmosphere. In some polynya regions, turbulent heat loss is as high as 300 W/m 2 such as along the Weddell Sea coast (see Morales Maqueda el a., 2004). However, mean turbulent heat loss in the study region was about 48 W/m 2 in 2018, which is much smaller than that in the St. Lawrence Island polynya (412 W/m 2 ; Pease, 1987) and Okhotsk Sea coast (471 W/m 2 ; Alfultis and Martin, 1987), 425 but comparable to the Northeast Water (NEW) polynya (31 W/m 2 ; Morales Maqueda el al., 2004). Nevertheless, new sea ice formation was somewhat slowed down in the study region because the air-sea temperature difference was much smaller than in the other regions listed above due to the anomalously warm air carried by southerly winds over the northern Greenland. 430 Finally, unlike the events in 2011 and 2017, sea ice in the northern Greenland region was not fully restored even a month after the peak ice removal on February 25 th , 2018 (Fig. 4d). Based on the RASM simulation between February 1 st to March 31 st , 2018, the study region lost approximately 208 km 3 of sea ice due to mechanical ice removal. Although sea ice was added later by the dynamic replenishment (102 km 3 ) and the thermodynamic ice growth (64 km 3 ), the deficit of SIV was about 42 km 3 , which is 435 equivalent to the 36 cm reduction of mean SIT over the study region. The CryoSat-2/SMOS data also confirmed the similar sea ice loss in the region; at the end of March 2018, the mean sea ice north of Greenland was relatively thinner (2.45 m), compared to the beginning of February (2.73 m). Hence, one could hypothesize that the 2018 winter polynya event could have contributed to the preconditioning of the polynya event in the following summer (Schweiger et al., 2021), which was observed at a similar 440 location. Hence the winter event might have been unprecedented by yet another measure, as polynya events have never repeated within a year over the study region except 2018.

Summary
Following the previous study of the 2018 winter polynya by Moore et al. (2018), this study demonstrates that all three observed winter open water events occurred under a specific wind pattern, 445 which removes sea ice mechanically out of the region. In February 2011 and 2017, the wind direction was primarily from the south and southwest (Figs. 4i and 4h). The size of polynya depended on the strength and persistence of southerly winds, and the polynya closure was a result of the relaxation of these conditions. The polynya was larger in February 2011 compared to the one in 2017 because wind duration was longer: 4 days of strong southerly wind (> 10 m/s) in 2011 compared to 2 days in 2017 450 (Fig. 8a). Hence, 31% more sea ice was removed in the former, although wind speed and wind stress were similar between the two years ( Table 1). The observed polynya in February 2018 was the largest one over the satellite SIC record since 1979. It resulted from much stronger and more persistent southerly-southeasterly winds (Fig. 4g) off the coast of northern Greenland, the direction which could also be more favorable for polynya opening (Kawaguchi et al., 2010). Although sea ice was 455 significantly thicker in the 1980s, more polynyas were not prevalent with thinner ice in the winters of 2015-2025 from the two RASM-DPLE ensembles. Given the similar rate of winter polynya occurrence and their apparent lack of sensitivity to the initial sea ice thickness, we conclude that a dominant cause of these winter polynyas stems from internal variability of atmospheric forcing rather than from the forced response to warming climate. 460 Table 1. Characteristics of the polynya events, when daily dynamic sea ice removal is greater than 10 km 3 /day (i.e., dynamic volume tendency (DVT) ≤ -10 km 3 /day) for more than three consecutive days, and the corresponding wind condition as well as sea ice thickness (SIT) over northern Greenland from the RASM hindcast simulation