Southern Ocean polynyas in CMIP6 models

Polynyas facilitate air-sea fluxes, impacting climate-relevant properties such as sea ice formation and deep water production. Despite their importance, polynyas have been poorly represented in past generations of climate models. Here we present a method to track the presence, frequency and spatial distribution of polynyas in the Southern Ocean in 27 models participating in the Climate Model Intercomparison Project phase 6 (CMIP6) and two satellite based sea ice products. Only half of the 27 models form open water polynyas (OWP ::::: OWPs), and most underestimate their area. As in satellite observations, 5 three models show episodes of high OWP activity separated by decades of no OWPs :::: OWP, while other models unrealistically create OWPs nearly every year. The :: In ::::::: contrast, ::: the : coastal polynya area in contrast is often overestimated : is :::::::::::: overestimated :: in :::: most ::::::: models, with the least accurate representations occurring in the models with the coarsest horizontal resolution. We show that the presence or absence of OWPs are : is : linked to changes in the regional hydrography, specifically the linkages between polynya activity with deep water convection and/or the shoaling of the upper water column thermocline. Models with 10 an accurate Antarctic Circumpolar Current (ACC) transport and wind stress curl have too frequent OWPs. Biases in polynya representation continue to exist in climate models, which has an impact on the regional ocean circulation and ventilation that require to :::::: should be addressed. However, emerging iceberg discharge schemes, vertical discretisation :::: more :::::::: adequate ::::::: vertical ::: grid :::: type : or overflow parameterisation, are anticipated to improve polynya representations and associated climate prediction in the future. 15

. The 27 CMIP6 models used in this study. Columns show the model names, ocean component, nominal horizontal resolution of the ocean (Ro [ ::: km]), atmosphere component, horizontal resolution of the atmosphere (Ra [ ::: km]), sea ice component, vertical discretisation scheme including number of vertical levels (z: depth-level, ρ: isopycnal, σ: terrain-following, several symbols: hybrid grid), data availability of monthly sea ice concentration (Cm), daily sea ice concentration (C d ), monthly sea ice thickness (Tm) and daily sea ice floe thickness (F T d  . :::::: Despite ::: the :::::: crucial ::: role ::: of ::::: OWPs : for sea ice production and deep water properties, an evaluation and comparison of their rep-60 resentation across current climate models is missing. Here we determine the characteristics, causes and impacts of Southern Ocean polynyas in 27 models that participated in the Climate Model Intercomparison Project phase 6 (CMIP6, Eyring et al., 2016). We present a new method to detect polynyas from daily and monthly sea ice concentration or thickness in CMIP6 models and observational data sets in Section 3. We use this approach to assess Southern Ocean polynya characteristics (spatial distribution, frequency and seasonality) for the 65 entire CMIP6 historical runand compare it : , ::::: check ::::::: whether :::::: major ::::: issues ::::::::::: documented :: in :::::: earlier ::::: CMIP ::::::: versions :::: still ::::::: remain, ::: and ::::::: compare ::: the :::::::: modelled :::::::: polynyas : with observational satellite data when available in Section 4. As observational data from within open water polynyas :::::: OWPs are very sparse, we continue with a comparison of the models' hydrography in sea ice covered conditions to that occurring during OWP events in order to assess OWPs' impact on the entire water column. We also analyse indicators of deep convection and upwelling as possible causes for the polynya formation and maintenance in CMIP6 70 models. Finally, with the obtained polynya statistics from the CMIP6 models, we determine whether the observed correlations between the SAM, Southern Hemisphere westerlies, ACC transport and the presence or magnitude of OWPs hold true across CMIP6 models (Section 5). We first consider the entire Southern Ocean south of 55°S, and then focus on the Weddell Sea region, between 65°W and 30°E.
2 CMIP6 output fields and observational data 75 We start with a description of the analysed data and then present the algorithm that we use in Section 3 to find polynyas.
We used for this study the daily sea ice concentration from OSI-450 , the daily sea ice SMOS thin ice thickness (Huntemann et al., 2014) :: two ::::::: remote ::::::::: observation ::::: based ::: sea ::: ice :::::::: products ::: for ::: sea :: ice :::::::::::: concentration ::: and :::::::: thickness :::: that :: we :::::::: introduce ::: in ::::::: Sections ::: 2.2 ::: and :::: 2.3, and the sea ice output from the historical run of 27 CMIP6 models (Eyring et al., 2016, listed in Table 1). We created this 80 2.1 ::::: Choice ::: of :::::: CMIP6 ::::::: models ::: We :::::: created : a : subset of all available CMIP6 models by filtering for models that participate in the 'historical' CMIP6 scenario, have at least one sea ice variable downloadable in a monthly or daily format and provide adequate projection and cell area information. The historical CMIP6 run covers the years from 1st January 1850 to 31st December 2014 and was : is : forced with observed historical greenhouse gas concentrations. As of the latest date of download (October 2020), only one ensemble 85 member was available for the majority of models. Consequently, we chose one representative ensemble member (r1i1p1f2 for CNRM models, r1i1p1f3 for HadGEM, r1i1p1f1 otherwise) for each model. This way, we ensure consistency with other CMIP sea ice area evaluations (Roach et al., 2020;Turner et al., 2013). For the assessment of polynya activity in the Southern Ocean, we use sea ice concentration and sea ice thickness data. In the 27 analysed models, the nominal horizontal resolution varies between 0.25°and 1.5°for the ocean (including the sea ice) and between 1°and 2.5°for the atmosphere (Table 1). 90 All computations were performed on the models' native grid unless specified otherwise. The grid cell area ('areacello') was used to compute the surface area of sea ice and polynyas. The models were purposely not detrended, as we want to determine 4 the accuracy of their historical run with ongoing climate change incorporated into the analysis. For the assessment of polynya activity in the Southern Ocean, we use the sea ice concentration and sea ice thickness as discussed below.

Sea Ice Thickness
For the CMIP models, two output variables related to the sea ice thickness are available. The sea ice volume per grid cell area ('sivol') has been available in earlier CMIP versions, and is commonly and hereafter referred to as equivalent sea ice thickness or just thickness. With CMIP6, a new sea ice variable called ice floe thickness ('sithick') is introduced and available for the majority of models (Table 1). The equivalent sea ice thickness sivol is available for 23 of the models, but only at a 120 monthly averaged resolution. Polynyas are dynamic processes, and can change their area and position drastically within few days (Nakata et al., 2015;Kern et al., 2007); thus it would be : is : optimal to analyse daily data from all data sources :: if ::::::: available.
The ice floe thickness sithick is available on a daily resolution for 21 of our models, but we found it not good ::: poor : for our polynya detection, as we show in Section 3.3.
For observational comparison of the sea ice thickness, we use the ::::: daily, satellite-based SMOS :::: "Soil ::::::: Moisture :::: and :::::: Ocean 125 ::::::: Salinity" ::::::: (SMOS) : thin sea ice thickness product (Huntemann et al., 2014)at a daily resolution. The product consists of maps of sea ice thickness up to 0.5 m, derived from the satellite-borne L-band radiometer SMOS. The data are distributed at a horizontal resolution of 12.5 km x 12.5 km and available from 1st June 2010 onwards. The sea ice thickness calculation from satellite data is a relatively new method that works best for ice thickness values of up to 50 cm, with an average retrieval error of about 30% of the retrieved value (Huntemann et al., 2014). It does not take differences in sea ice concentration into account and might thus 130 be biased in regions with low sea ice concentration. The accuracy of the method is negatively affected by melt ponds (ponds that the form on the sea ice in the melting season), but low air humidity causes melt ponds to occur much less frequently on the Antarctic sea ice compared to the Arctic sea ice (Andreas and Ackley, 1982). The thin sea ice retrieval is a relatively new method that is more established for Arctic regions (e.g. Tietsche et al., 2018), but papers about its quality and applicability for the Southern Ocean are in preparation (Mchedlishvili et al., 2021;Tian-Kunze et al., 2021) ::::::::::::::::::::::::: (e.g. Mchedlishvili et al., 2021) 135 . Even though the 10 years time period is too short for deriving climatological mean values, we decided to include it as an observational baseline, where relevant, to compare the CMIP6 data to.
Due to the described advantages and disadvantages of using either sea ice concentration (siconc) or thickness (sivol), we have decided, where practical, to analyse both and present the resulting data side by side in sections 3 and 4.

140
We do not limit our analysis to the surface of the ocean, but also analyse vertical ocean profiles within CMIP6 polynyas. The majority of models use a z-level vertical grid in the ocean, with the exceptions of GFDL-CM4/ESM4, MIROC6/ES2L and NorCPM1. The GFDL models use isopycnal coordinates in the interior ocean and rescaled ::::: rescale :: to : geopotential vertical coordinates in the mixed layer. The MIROC6 model uses hybrid z -sigma coordinates between the sea surface and a fixed geopotential depth and a z-level grid below. NorCPM1 uses isopycnal coordinates in the interior ocean and z coordinates in 145 the mixed layer; it is also the only model using data assimilation . To investigate the hydrography we extracted vertical profiles from the monthly ocean salinity ('so') and potential temperature ('thetao'). For observational comparison of the water properties, data from a SOCCOM vertical profiling float (Johnson et al., 2018) is used. We use temperature and salinity profiles measured by the SOCCOM float 5904471, because it provides a long record of profiles including rare vertical profiles from within the 2017 Maud Rise OWP .

Averaging methods
For readability, we condensed the data further in several ways depending on the objective of the analysis. To analyse the spatial distribution of polynyas, we present sea ice maps for September with polynya occurrences highlighted ( Fig. 3, A3, A4). We 225 chose the month of September as it typically features the austral sea ice maximum and is therefore discussed in related CMIP6 sea ice and Southern Ocean evaluations (e.g., Beadling et al., 2020;Roach et al., 2020). However, we also give the different seasonality of the modelled coastal and open water polynyas ::::: OWPs in section 4.2. In Figures 6, 7, A1, A2 and for the polynya areas in Table 2, we used equation 1 first and took the mean over all model years then. For the sea ice extent in Table 2 for daily and monthly variables respectively. For the figures within this paper we present the result of the analysis of either the sea ice concentration or thickness in either daily or monthly resolution to keep the data comparable. Not all plots include all 27 analysed models, since some models do not provide all output variables (Table 1). To avoid unnecessary repetitions in the paper, we do not always present the results from all sea ice variables used, but keep some results exclusive to the supplementary 235 figures.
In Section 4, we limit the analysis to the period May to November in order to filter out summer polynyas. Throughout the paper, we use the mean yearly polynya area for the winter season if not specified otherwise. In sea ice melting summer conditions, polynyas become often very large, but we want to focus on winter polynyas due to their importance in heat exchange and ice and deep water formation. Moreover, it is not recommended to use the SMOS thin ice product during melting season 240 (Huntemann et al., 2014). When comparing model output data to the observational sea ice concentrations, we further limit the yearly mean values to a common time period from 1979-2015 to ensure consistency.
A combined approach in which either low sea ice concentration or thickness or the product of both classifies an area as polynya would be possible, but out of the 21 models that include daily sea ice thickness, three were lacking daily sea ice con-

Caveat
To asses and compare polynya activity in a large number of models, we analyzed different output variables for equivalent sea ice thickness and concentration. The equivalent sea ice thickness 'sivol' was only available at a monthly resolution. We hope to see more daily sea ice data for all models published in the next CMIP iteration.

265
The current algorithm differentiates coastal and open ocean polynyas by their direct connection to the continent. In general, we find a good agreement between the algorithm's classification and our visual validation. However, occasionally the algorithm classifies a polynya that is (at least partially) located on the coastal shelf as an open water polynya; or a large open water polynya neighbouring land on one grid cell, as a coastal polynya. While this classification holds true for a strict definition of the polynyas, a polynya on the coastal shelf is unlikely to be a sensible heat polynya driven by deep water convection.

Polynya statistics in CMIP6
300 The aim of our study is to determine whether the representation of polynyas in CMIP6 models is accurate in terms of location (section 4.1) and frequency (section 4.2). To do so, we evaluate ::::::: compare the observational products and the model output with the same methods introduced in Section 3. We finally investigate the effect of the modelled polynyas on modelled stratification in section 4.3.
In agreement with observations, the models show small coastal polynyas during winter that grow in size and in time, merging 315 with other polynyas in spring/summer (Nov-Jan, see  (Table 1). Seven 320 of these families show polynyas at very similar locations for all members inside one family (Fig. 3, A3, A4); the only exception is the ACCESS models. This indicates that the position of polynyas is mostly determined by the model properties and less by coincidence of sporadic polynya occurrence or long term variability. Comparison of the results from sea ice thickness (Fig. 3) and sea ice concentration (Fig. A4, A3) shows similar locations, but we found on average almost twice the polynya area by :::: with the sea ice thickness threshold method (Table 2).  Most models (15 when using thickness, 21 concentration) show a larger fraction of coastal polynya than OWP area, while in the rest the OWPs prevail as a larger fraction (Fig. 6). Surprisingly ::::::::::: Surprisingly, the model ACCESS-CM2 shows only 345 coastal polynyas while ACCESS-ESM1-5 shows an overabundance of OWPs. These two models share similar versions of the ACCESS-OM2 ocean model component and CICE sea ice model component (Table 1). We further investigate the cause for this difference and its effect on the water stratification in section 5.3.

Polynya frequency and seasonality
In Section 4.1, we discussed the average polynya areas over 35 years. A closer look at individual yearly values from the full 350 historical run reveals that the polynya areas vary from year to year (Figures 7); for the coastal polynyas the ratio between the year with the largest polynya area to that with the smallest one is at least 2.5 (lowest for MRI-ESM2, highest for CESM2-FV2, HadGEM3-GC31-LL, MPI-ESM-1-2-HAM, SAM0-UNICON, which have one ore more years without any coastal polynyas).
The distribution of the OWPs has one fundamental difference to the one for the coastal polynyas. All models except for the 365 BCCs show years without any or with negligible OWP activity (<1000 km :: 10 ::: 000 :::: km). The OWP activity we have found in Section 4.1 for the EC-Earth3, GFDL and MPI models is mostly caused by some isolated, large scale OWP events while the majority of years have little to no OWPs. OWPs in the models are not randomly distributed in time but instead regularly re-open for several consecutive years, as seen in the MPI, ACCESS, BCC, CAMS-CSM1 and EC-Earth models (Fig. 8, A6, A7). This is interesting since the sea ice in the Weddell Sea melts almost completely every summer. Thus polynya formation must either For the majority of models, the coastal polynyas appear randomly distributed over the years (Fig. 8, A6, A7). In agreement with observations, the models show the largest coastal polynyas between November and January (Fig. 9b,c). However, the models underestimate the coastal polynya area in November and December followed by an equally strong overestimation 380 for January and February (Fig. 9b). The season with the highest OWP activity in the CMIP6 models is between August and November (Fig. 9c,f), while the the sea ice observations show a local maximum in polynya occurrence during ice formation in June and an absolute maximum during ice melt in November/December. The seasonality of the OWPs does not agree well between CMIP6 and observations.
If the stratification was already weak or the density at the surface is increased by brine rejection of newly formed sea ice, deep convection may commence (weakly stratified cold profiles, Fig. 10).
Re-arranging the profiles into Temperature-Salinity diagrams (right of Fig. 10 ::::: d,h,l,p), the core of the comparatively warm Circumpolar Deep Water (CDW), at approximately 1000 m depth, is easily visible from the temperature maximum. Under 450 sea ice, the water masses on either side of the CDW are clearly separated: the surface layer is cold and fresh, while the deep water is saltier and warmer than the surface. The pycnocline between these layers is at a depth of 100 to 200 m, except for the BCC models, which show very weak stratification and no clear water mass layers can be identified. Inside the polynyas, the differentiation of the water masses is almost gone (Fig. 10). The models with the highest OWP activity (ACCESS-ESM1.5, BCC and MPI) have low ::::: values ::: of ::: the :::::::::::: Brunt-Väisälä :::::::: frequency ::::::: squared : (N²values : ) : even under sea ice, where they hardly 455 exceed 0.5 10 −5 s −2 (Fig 10 ::::: c,g,k,o). The deep mixing is transporting :::::::: transports : heat and salt into the surface layer, reducing its density difference compared to deeper layers. This weakened stratification can remain in the water for a significant period of time, preconditioning the region for new polynyas to form in following years (Martinson et al., 1981).

Discussion
In the previous sections, we looked at the representation of polynyas in CMIP6 models and their effect on the local water 460 stratification. We found large variations among models, with the majority underestimating OWP area but a few largely overestimating it. The hot spot :::::: hotspot for OWPs in CMIP6 and observations is the Weddell Sea. Now we want to turn our attention to finding the causes for the difference in polynya activity in the CMIP6 models. We start with a comparison of Climate models versus Earth System models and then discuss whether the model's winds and Antarctic Circumpolar Currents can explain its OWP representation. Then we will examine the possible reasons for the ::::: strong overestimation of OWPs in many :::: some : models 465 and give an outlook about how we expect these biases to change in the near future.
Open water polynya activity is commonly attributed to the Southern Annular Mode ::::: SAM, the strength of the Southern Hemisphere winds and the transport of the ACC (e.g. Swart et al., 2018;Cheon et al., 2014;Campbell et al., 2019;Behrens et al., 490 2016). We combine our polynya characteristics from section 4 with the statistical evaluation of Southern Ocean properties from Beadling et al. (2020) to determine potential across-model relationships between Southern Ocean properties and polynya activity. Since Beadling et al. (2020) presents data from 1986-2005, we restrict our analysis to the same period and the 23 models we have in common. The CMIP6 across-model correlation between the total wind stress forcing of the westerlies (Beadling et al., 2020) with OWP activity (Table 2) is 0.4 ::: 0.49: that is, the models with the strongest Southern Hemisphere
We have shown that overabundant OWPs are a common issue in models (sections 4.1 and 4.2); that is a problem because it affects the whole water column (section 4.3). We have now found relationships with the models' wider representation of the 530 Southern Ocean, but these relationships have some exceptions. In a case study, we will now investigate the specific case of two seemingly similar models that have extremely different OWP representations, ACCESS-CM2 and ACCESS-ESM1.5.

Overestimation of OWPs in models and future perspectives
We found that most models overestimate polynya activity in comparison to the observations (Table 2). We believe this to be general known ::::::::: knowledge in the community given that within model documentation, there are many descriptions of approaches 535 to reduce deep convection and the presence of OWPs in the Southern Ocean (listed and highlighted with a '*' in supplementary   Table A1). We first describe approaches to reduce OWPs and conclude with future perspectives.

545
One of the development foci for the GFDL models was to reduce Southern Ocean polynyas . For the GFDL-CM2 models, the formation of super polynyas in the Southern Ocean was addressed by increasing the near-infrared albedo of glaciers and snow covered ice caps in order to increase coastal freshwater inflow. This is reported to delay the formation of super polynyas in the model, but not to prevent them completely . ::: One :::::: reason ::: for ::: this :::::: partial The CESM2 models feature an overflow parameterisation (Briegleb et al., 2010), which can transport dense bottom water down the shelf and should help with more realistic bottom water formation. This can lead to a more stable stratification, but Heuzé (2015) found that it also caused an overestimation of coastal polynyas in their CMIP5 model version. In our study, we 555 find that the CMIP6 CESM2 models show coastal polynyas, but never OWPs. The observed polynya area is underestimated by at least 30% by all analyzed CESM2 models.
We found that both ACCESS models show more than double the polynya area of the observational data (Fig. 6, Fig. 9).
Even though the ACCESS-ESM1.5 and the ACCESS-CM2 models share most model components, the former shows prevalent OWPs while the latter has mostly coastal polynyas. In the assessment of ACCESS-ESM1.5, Bi et al. (2013) found that OWPs suppresses the OWP activity in ACCESS-CM2 (Fig. 3). The freshwater transport of icebergs somewhat suppresses too frequent open ocean deep convection events (Heuzé, 2021), but ACCESS-CM2 forms too large coastal polynyas instead, much 565 more frequent than ACCESS-ESM1.5. The extent of these coastal polynyas is unrealistically high and the problem of overall too large polynya areas remains. The suppressed OWP activity due to freshwater discharge from icebergs affects the whole water column in the Weddell Sea (Fig. A8). In ACCESS-CM2, the Weddell Deep Water is 0.7°C warmer and the maximum temperature of the CDW is 1.1°C warmer than in ACCESS-ESM1.5.
Here we discussed some example :::::::: examples of modellers' approaches to reduce deep convection and OWPs in the Southern

570
Ocean. In the model documentations (listed in supplementary Table A1), we found more examples with the same aim. On the other hand, several models have found a way to prevent this :::: these : issues. The CESM2 models have an overflow parameterisation (Briegleb et al., 2010) and do not show any OWPs. This parameterisation can provide an effective way for realistic bottom water formation, until higher model resolutions and further model improvements (e.g. better vertical discretisation schemes) allow for actual bottom water formation in coastal polynyas. Also, isopycnal coordinates are beneficial for the representation of down 575 slope flows, the deep water formation in coastal polynyas. This in turn results in more realistic deep water properties (Heuzé et al., 2013). We found that no model with isopycnal coordinates (Table 1) shows an unrealistic amount of OWPs (Table 2 and Ocean deep water properties and OWPs when these models successively upgrade to MOM6.
In section 4.3 we ::: We have shown how OWPs significantly reduce the local stratification. We found indications of upwelling as well as deep convection, which eventually change the long term water properties in the Weddell Sea, notably removing the stratification and hence making future polynyas more likely. In section :::::: Section 5.3 we discussed the long term effect of coastal versus open water polynyas in the Weddell Sea using two different ACCESS model versions.

32
Appendix A: Supplementary Figures and Tables Figure A1. Correlation plots for the polynya areas computed from daily (x) vs monthly (y) sea ice concentrations as presented in Table 2.
Each dot represents the averaged polynya area of one year; red are the coastal polynyas, blue the open water polynyas ::::: OWPs : and black the sum of both. Where applicable (e.g. polynyas present, p-value < 0.05) the legend includes the slopes of the correlation in the order coastal/open water/combined polynyas. Because the area of polynyas varies by more than one order of magnitude, we normalized the axes of each plot by the highest numerical value for each model and give its value in the lower right corner of each plot. Figure A2. Correlation plots for the polynya areas computed from monthly sea ice thickness (x) vs monthly sea ice concentration (y) as presented in Table 2. Each dot represents the averaged polynya area of one year; red are the coastal polynyas, blue the open water polynyas ::::: OWPs and black the sum of both. Where applicable (e.g. polynyas present, p-value < 0.05) the legend includes the slopes of the correlation in the order coastal/open water/combined polynyas. Because the area of polynyas varies by more than one order of magnitude, we normalized the axes of each plot by the highest numerical value for each model and give its value in the lower right corner of each plot.  36 Figure A5. The seasonal cycle of sea ice concentration (blue colours) and polynya probability (red-yellow) for the GFDL-CM4 model.
The 28 CMIP6 models used in this study; full name and references with model description.
The references marked with a "*" contain information about deep convection in the Southern Ocean. N/A indicates that no paper has been published yet for the CMIP6 configuration. Model Name Descriptive Name