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 (OWPs), and most underestimate their area. As in satellite observations, three models show episodes of high OWP activity separated by decades of no OWP, while other models unrealistically create OWPs nearly every year. In contrast, the coastal polynya area 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 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 an accurate Antarctic Circumpolar Current 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 should be addressed. However, emerging iceberg discharge schemes, more adequate vertical grid type or overflow parameterisation are anticipated to improve polynya representations and associated climate prediction in the future.

1 Introduction

Polynyas are areas of open water surrounded by sea ice. They are common features within the Southern Ocean winter sea ice and are often classified into two different categories: coastal polynyas and open-water polynyas (OWPs). Coastal polynyas are usually latent heat polynyas that are kept open by winds that drive the sea ice away from the coastline or an obstacle, like an iceberg (Morales Maqueda et al., 2004). OWPs are kept open by thermodynamic processes, such as upwelling of sensible heat or diffusive fluxes (Martinson and Iannuzzi, 1998) that transport warm water masses upwards, melting the sea ice and keeping the regional ocean ice-free. The heat source that keeps OWPs from freezing over is the presence of comparatively warm and salty Circumpolar Deep Water (CDW), which is usually located just below the base of the upper mixed layer (Santoso et al., 2006). These warmer waters transport heat from the depth to the surface by free or forced convection, caused by surface cooling, brine rejection, or wind and shear stresses respectively (Williams et al., 2007). While open-ocean deep convection is closely linked to OWPs (Cheon and Gordon, 2019), deep convection is not a sufficient condition for OWPs to occur (Dufour et al., 2017).

Polynyas play a key role in sea ice and deep water formation. In winter, the strong temperature contrast between the open water in the polynyas and the cold Antarctic air causes significant heat loss of several hundred W m⁻² from the ocean to the atmosphere (Willmott et al., 2007). In fact, while the area of coastal polynyas is only about 1 % of the sea ice area in the Southern Ocean, 10 % of Antarctic sea ice is produced there owing to this intense cooling (Tamura et al., 2008). The rapid ice production in coastal polynyas also produces a water mass of extremely high density, the Dense Shelf Water. This water is considered to be the main source of Antarctic Bottom Water, which forms the deepest layer of the global oceans (Orsi et al., 1999).

Table 1 The 27 CMIP6 models used in this study. Columns show the model names, ocean component, nominal horizontal resolution of the ocean (R_o [km]), atmosphere component, horizontal resolution of the atmosphere (R_a [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 (C_m), daily sea ice concentration (C_d), monthly sea ice thickness (T_m) and daily sea ice floe thickness (FT_d).

Some coastal polynyas open at the same location every year (e.g. the Ross Sea polynya), while OWPs are observed only once per decade or less. The Weddell Sea polynya has been the largest OWP observed to date and affected the water properties in the Weddell Sea for decades (Zanowski et al., 2015). It was first observed in the winters of 1974–1976, with an average area of up to 300 000 km² (Carsey, 1980); produced about 5 Sv of Antarctic Bottom Water (Wang et al, 2017); and led to a significant cooling at mid-depth (Cheon and Gordon, 2019). Only minor polynyas have been observed in the region from 1976 to 2016 (Cheon and Gordon, 2019), when the Maud Rise polynya (Swart et al., 2018; Heuzé et al., 2019; Francis et al., 2019) reached an area larger than 50 000 km², similar to the polynya of October–November 1973 (Cheon and Gordon, 2019). These Weddell Sea and Maud Rise polynyas were preceded by a phase of positive Southern Annular Mode (SAM) anomalies and strong Southern Hemisphere westerlies (Gordon et al., 2007; Campbell et al., 2019; Cheon et al., 2014). Behrens et al. (2016) found that deep convection, as it happens in OWPs, is strengthening the meridional density gradients and increases the Antarctic Circumpolar Current (ACC) transport on multidecadal timescales.

Polynya detection in observational data (e.g. Markus and Burns, 1995; Kern et al., 2007; Ohshima et al., 2016) and ocean reanalysis products (Aguiar et al., 2017) and the role of OWPs in spurious deep water formation in CMIP5 (Heuzé, 2015) have been discussed before. The formation of OWPs in models has been shown to be very sensitive to vertical mixing parameters (e.g. Kjellsson et al., 2015; Heuzé, 2015), initial sea ice conditions (e.g. Kjellsson et al., 2015), stratification (e.g. De Lavergne et al., 2014; Stoessel et al., 2015) and model resolution (e.g. Dufour et al., 2017; Kurtakoti et al., 2018; Lockwood et al., 2021). A weak background stability, especially in the Weddell Sea (Wilson et al., 2019), makes the ocean susceptible to convective overturning due to model inaccuracies such as a lack of dense shelf water overflows (Dufour et al., 2017) or heat buildup due to insufficient vertical mixing (Heuzé et al., 2015). Spurious OWP appearance or deep convection in the Weddell Sea remains a challenge in modern climate models (e.g. Held et al., 2019; Sellar et al., 2019; Mauritsen et al., 2019; and references marked with an asterisk (*) in Table A1). Ocean convection due to static instabilities is an important process in the formation of OWPs, which is not modelled directly due to the relatively coarse resolution of many CMIP6 models (Table 1) but parameterised instead (e.g. Hasumi, 2000; Madec et al., 2017). Despite the crucial role of OWPs for sea ice production and deep water properties, an evaluation and comparison of their representation 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 Sect. 3. We use this approach to assess Southern Ocean polynya characteristics (spatial distribution, frequency and seasonality) for the entire CMIP6 historical run, check whether major issues documented in earlier CMIP versions still remain and compare the modelled polynyas with observational satellite data when available in Sect. 4. As observational data from within 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. 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 (Sect. 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

We start with a description of the analysed data and then present the algorithm that we use in Sect. 3 to find polynyas. We used for this study two remote observation-based sea ice products for sea ice concentration and thickness that we introduce in Sects. 2.2 and 2.3, as well as the sea ice output from the historical run of 27 CMIP6 models (Eyring et al., 2016, listed in Table 1).

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 1 January 1850 to 31 December 2014 and is forced with observed historical greenhouse gas concentrations. As of the latest date of download (October 2020), only one ensemble member was available for the majority of models. Consequently, we chose one representative ensemble member (r1i1p1f2 for CNRM models, r1i1p1f3 for HadGEM and 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). 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). 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 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.

2.2 Sea ice concentration

Sea ice concentration (CMIP6 parameter: “siconc”) is available at a daily resolution for only 18 of the 27 models, whereas it is available at a monthly resolution for all 27 models (Table 1). Furthermore, the daily sea ice concentration is not available for the full historical period for the models MRI-ESM2-0 (1 January 1920–31 December 2015) and SAM0-UNICON (1 January 1950–31 December 2015). Routine satellite-based sea ice concentration observations are available from January 1979 onwards, so there are 35 years of overlap with the historical CMIP6 model run on which we can perform our comparisons. We use the daily satellite-observation-based Global Sea Ice Concentration climate data record (OSI-450) (Lavergne et al., 2019) for the time period 1979–2015. For the time after 2015, an extension with the name OSI-430-b is available that is processed with the same algorithms (Lavergne et al., 2019). OSI-450/430b are sea ice concentration products, computed from SMMR (1979–1987), SSM/I (1987–2008), and SSMIS (2006–2020) microwave radiometers and ECMWF ERA-Interim reanalysis data (Dee et al., 2011). They are provided at a 25 km × 25 km horizontal resolution and daily temporal resolution from 1 January 1979 to 2020 onwards (every other day until 1985). The uncertainty in sea ice concentration is less than 4 % on average (Lavergne et al., 2019).

Using sea ice concentration to detect polynya activity has the advantage that it is available for the largest number of CMIP6 models (Table 1) and has been observed by satellite for more than 40 years. The disadvantage is that sea ice concentration is a poor choice to detect coastal polynyas that are often covered by newly forming thin sea ice in winter, and their full extent is better characterised by a low sea ice thickness (Ohshima et al., 2016). Polynya detection by sea ice thickness thresholds is well established (e.g. Ohshima et al., 2016; Nakata et al., 2015; Kern et al., 2007) and detects even polynyas that are covered by a thin layer of newly formed sea ice. 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 Sects. 3 and 4.

2.3 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 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 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 poor for our polynya detection, as we show in Sect. 3.3.

For observational comparison of the sea ice thickness, we use the daily, satellite-based Soil Moisture and Ocean Salinity (SMOS) thin sea ice thickness product (Huntemann et al., 2014). 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 × 12.5 km and available from 1 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 be biased in regions with low sea ice concentration. The accuracy of the method is negatively affected by melt ponds (ponds that 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 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 (e.g. Mchedlishvili et al., 2021). Even though the 10-year 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.

2.4 Vertical ocean profiles

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 rescale to geopotential vertical coordinates in the mixed layer. The MIROC6 model uses hybrid z–σ 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 the mixed layer; it is also the only model using data assimilation (Counillon et al., 2016). 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) are 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 (Campbell et al., 2019).

3 Polynya detection

In this work, we present a new method to detect and distinguish coastal and OWPs using either sea ice concentration or equivalent sea ice thickness (Fig. 1). We apply this method to the observational, satellite-based products and the CMIP6 models listed in Table 1.

Figure 1 A four-step algorithm classifies surface type areas in the marginal ice zone. (a) August 2002 monthly mean sea ice concentration from ACCESS-ESM1.5. (b) Sea ice map after covering the open ocean with the flood fill algorithm. (c) Application of a threshold filter, highlighting the total polynya areas. (d) Flood fill algorithm is applied to Antarctic fast land; the remaining open water or low sea ice areas are classified OWPs. The ice shelves appear in light grey and the ocean in dark grey.

Figure 2 Illustration of the effect of different sea ice concentration thresholds on the detection of polynyas. (a) Observational sea ice concentration (OSI450) from 25 November 2017. (b) Combined area of coastal polynyas (orange) and OWPs (green) for sea ice concentration thresholds from 0 % to 100 %. The vertical lines correspond to the results of our algorithm for (c) 15 %, (d) 30 % and (e) 60 % sea ice threshold.

3.1 Our algorithm to detect polynyas

With the aim of detecting polynyas, we start with the sea ice concentration or thickness (Fig. 1a). To mask out the open ocean beyond the northern sea ice extent, we use a “flood fill” algorithm from the scikit-image library (Van der Walt et al., 2014). Starting from a grid cell with no sea ice, the seed, the algorithm detects similar cells below a specified sea ice concentration/thickness and masks them out, effectively “filling them” with ice (Fig. 1b). Afterwards, a maximum sea ice threshold filter returns all grid cells that are classified as “polynyas” (Fig. 1c). To differentiate coastal polynyas from OWPs, in a third step the flood fill algorithm is applied to the Antarctic continent (Fig. 1d). The Antarctic continent and ice shelves are represented by zero sea ice concentration/thickness in most models; for the exceptions that mark these grid cells in another way (e.g. infinite or negative value), we map this values to zero first. Coastal polynyas, adjacent to the land or ice shelves, are hence also covered by the flood fill algorithm. The grid cells below the threshold that remain are classified as “OWP”. The polynya area is the sum of the grid cell areas. The obtained OWP areas can be subtracted from the total polynya area to obtain the coastal polynya areas.

Depending on the chosen threshold and sea ice variable, the resulting polynya and sea ice areas differ. For example a closed layer of thin sea ice will be detected as a polynya if using the threshold thickness criteria but not if using the threshold sea ice concentration (siconc) criteria; an area with relatively sparse but thick sea ice will be classified as a polynya by the threshold sea ice concentration (siconc) criteria but not if using the threshold sea ice thickness (sivol) criteria. In the literature, a wide range of values for the concentration thresholds are in use. Previous studies on polynyas tend to choose a high sea ice concentration threshold to maximise the number of polynyas detected (e.g. Arbetter et al., 2004; Gordon et al., 2007), while sea-ice-specific studies choose a low threshold to capture all areas with sea ice. For example, Shu et al. (2015) and Beadling et al. (2020) use a 15 % threshold to estimate the sea ice extent in CMIP5 and CMIP6 models, Kern et al. (2007) classify areas with up to 45 % sea ice coverage as open-water areas, and Smedsrud (2005) points out that sea ice coverages of 80 %–90 % already allow for significant heat exchange and can therefore be classified as polynyas. The influence of different sea ice thresholds on our four-step algorithm is visualised in Fig. 2: the areas that are classified as polynyas increase with higher sea ice concentration thresholds, until the first polynya areas merge with the open ocean and thus become embayments (Fig. 2b, e). High sea ice concentration thresholds (up to ∼ 85 %) maximise the number of detectable coastal polynyas but lead to poor detection of OWPs (Fig. 2e). The higher the chosen sea ice threshold, the lower the area that is classified as sea ice. We found that almost no coastal polynyas are detected at a 15 % threshold in our observational product (Fig. 2b, c) and in most CMIP6 models (not shown). Therefore, we chose a 30 % sea ice concentration threshold as a good compromise: higher thresholds lead to a strong negative bias in sea ice area compared to other papers, and polynyas can become embayments as shown in Fig. 2e; lower values leave too many polynyas undetected. For thresholds higher than 60 %, the OWP area is increasing again (Fig. 2b) after the big OWP in the Weddell Sea merged with the open ocean. However, this signal is relatively noisy in time (not shown) and sensitive to small threshold changes (Fig. 2b). While more and more polynyas become embayments, new areas with more than 60 % sea ice concentration are classified as “polynyas”, and the area classified as “sea ice covered” in between shrinks. For the sea ice thickness threshold, we chose a value of 12 cm. This is consistent with Nakata et al. (2015), who give the total polynya extent as the areas occupied by open water and thin ice (up to 12 cm). For a discussion of the effect of different sea ice thickness thresholds, see Nakata et al. (2015) and Kern et al. (2007).

Some CMIP6 models provide only monthly averaged data, others include the sea ice concentration at daily resolution (Table 1). It is obvious that monthly data are unsuitable to detect a phenomenon that may last but a few days and that consequently the usage of monthly averaged sea ice output may affect the accuracy of the mean polynya areas. When comparing the daily observational data to monthly model output, we computed the monthly time mean of the sea ice concentration and thickness first and applied our algorithm then. To provide yearly averaged winter polynya area for each model and observational product, we used yearly polynya area averages:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& \stackrel{\mathrm{‾}}{A_y}\left(\text{daily data}\right)={\displaystyle \frac{\sum _{\text{1 May}}^{\text{15 Nov}}{A}_{\mathrm{d}}\left(t\right)}{N_{\mathrm{days}}}},\\ & \stackrel{\mathrm{‾}}{A_y}\left(\text{monthly data}\right)={\displaystyle \frac{\sum _{\text{May}}^{\text{Nov}}{A}_{\mathrm{m}}\left(t\right)}{N_{\mathrm{months}}}}.\end{array}\end{array}$$

In Eq. (1), A_d(t) and A_m(t) are the polynya areas derived from 1 d or 1 month. N_days and N_months are the total number of days or months in the time series. As we show in Sect. 4.2, coastal polynyas grow very large in the summer season, but we focus on winter polynyas here due to their large importance in ocean–atmosphere heat exchange and sea ice formation. We used the results from Eq. (1) to create scatter plots with least squares linear regressions (Figs. A1 and A2). For the observational sea ice concentration, the slope of the regression fit is 0.73, indicating that total polynya areas are underestimated by 27 % if computed from monthly averages and should be derived from daily data wherever possible. The regression is worst with OWP, probably because of the small number of OWPs in the observational data (Fig. A1). Seven of 18 CMIP6 models have a regression slope higher than or equal to the observational product (ACCESS-CM2, BCC-CSM2-MR, CNRM-CM6-1, CNRM-ESM2-1, EC-Earth3(-Veg), MPI-ESM1-2-HR), i.e. a better agreement between daily and monthly polynya areas, while the rest have a lower slope. We find a higher correlation for models including larger polynyas (total polynya area $>\mathrm{100}\times {\mathrm{10}}^{\mathrm{3}}$ km²) compared to models that have only minor polynyas. We suspect that smaller polynyas have a shorter average life span and are more prone to relative changes in area due to ice drift.

The correlation of yearly polynya areas from equivalent sea ice thickness (sivol) with sea ice concentration (siconc) shows a slope of 0.31 for the observational data sets SMOS thin ice and OSI-450/430 (Fig. A2). Here we could only use the common time period from 2010 to 2020. The low slope indicates that our algorithm classifies 69 % less area as polynya when using sea ice concentration. This is not necessarily an erratic underestimation in this case but rather originates in different definitions of a polynya either using a concentration or thickness threshold. However, for 12 of the 23 models this slope is higher, indicating a closer relationship between low sea ice thickness and concentration. This is likely due to the conceptual differences between the observational sea ice thickness and the CMIP6 sea ice thickness. The observational sea ice thickness retrieval algorithm assumes sea ice concentrations of 100 % (Huntemann et al., 2014), while the sea ice thickness variable we use for our analysis is an area-averaged sea ice thickness and is thus directly connected to changes in sea ice concentration. We found low slopes in models with more coastal polynyas (e.g. MRI-ESM2.0: 0.00, IPSL-CM6A: 0.05, CESM2: 0.11) but also in the BCC models that form wide areas of thin sea ice (<10 cm) followed by a thicker ring along the outer sea ice edge (not shown).

3.2 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 (Figs. 3, A3, A4). We 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 OWPs in Sect. 4.2. In Figs. 6, 7, A1, and A2 and for the polynya areas in Table 2, we used Eq. (1) first and took the mean over all model years then. For the sea ice extent in Table 2 and the maximum polynya extents in Fig. 6, we used

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& \stackrel{\mathrm{‾}}{A}={\displaystyle \frac{\sum _{\mathrm{years}}\max \left(A_{\mathrm{d}}\left(t\right)\right){|}_{\text{1 May}}^{\text{15 November}}}{N_{\mathrm{years}}}},\\ & \stackrel{\mathrm{‾}}{A}={\displaystyle \frac{\sum _{\mathrm{years}}\max \left(A_{\mathrm{m}}\left(t\right)\right){|}_{\mathrm{May}}^{\mathrm{November}}}{N_{\mathrm{years}}}}\end{array}\end{array}$$

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 Appendix figures.

In Sect. 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 often become 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 (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.

3.3 Sea ice floe thickness

The sea ice floe thickness (sithick) is defined as the actual thickness of the sea ice in the CMIP6 models, averaged over the ice-covered part of the grid cell (Huntemann et al., 2014). While the detected polynya locations and areas in our analysis of the sea ice concentration and thickness (Figs. 3, A3, A4) agree reasonably well with observations for most models and will be discussed further, we could not achieve satisfactory results if using only the sea ice floe thickness as input data (Fig. 4). The majority of the models (e.g. ACCESS, BCCs, CESM2s, MPIs, IPSL-CM6A-LR, MRI-ESM2-0, SAM0-UNICON) included increasing sea ice floe thickness beyond the outer sea ice extent of the concentration maps. These models also have a considerable mean floe thickness inside areas that would be classified as polynyas by low sea ice concentration and equivalent thickness (e.g. Fig. 4j, k). When using our thickness threshold algorithm with this data set, we obtain large OWPs in areas that hardly correspond to the observations (Fig. 4a–i), because the high ice floe thickness anomalies at the outer sea ice edges encircle large areas with thinner ice floe thickness.

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, 3 were lacking daily sea ice concentration. Where daily data are provided, they are sometimes only available for a limited period (MRI-ESM2-0, SAM0-UNICON). Moreover, we found that the MPI models show a strange behaviour, where the floe thickness is either > 0.5 m or zero but never takes values < 0.5 m. We assume that the differences in implementation of the floe thickness output variable are due to the novelty of this variable with CMIP6, so that no conventions over the exact interpretation and implementation of the variable have formed yet. We will focus instead on the equivalent sea ice thickness (sivol) and concentration (siconc) variables for the rest of the paper.

3.4 Caveat

To assess and compare polynya activity in a large number of models, we analysed 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.

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 it classifies 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. For example, according to our analysis, IPSL-CM6A-LR has few OWPs but mostly coastal polynyas; yet these coastal polynyas reach far out of the continental shelf, and this model was reported to undergo deep convection (Heuzé, 2021), which suggests a misclassification of some OWPs as coastal polynyas. Even though such polynyas adjacent to the coast are strictly speaking coastal polynyas and classified as such by our algorithm, the underlying physical process is likely deep water convection. An improvement of our method could be achieved by taking the water depth at each grid cell into account in the classification algorithm. Thus polynyas located on the shelf could be classified as coastal polynyas, while coastal polynyas extending further into the open ocean could be (partially) classified as OWPs.

A possible negative bias in our polynya area estimations could be sea ice embayments. Large embayments from the open ocean into the sea ice can reach latitudes and areas comparable to those of OWPs, especially in the melting season. We use a strict definition of polynyas as areas surrounded by sea ice and do not account or compensate the results for eventual embayment areas. Another bias could be the choice of one ensemble member for each model. In several models (GFDL-, ACCESS-, UKESM1-0-LL) multi-centennial variability causes periods with high polynya activity followed by times with few polynyas (Sellar et al., 2019). The ACCESS-ESM1-5 shows a long-term warming trend of more than 1 ^∘C over the 165 model years (not shown). The BCC models have a strong multidecadal variance in sea ice extent (Beadling et al., 2020), which also affects the polynya areas (Sect. 4.2). In these models with long-term variability, the results for the polynya activity could be less representative because we derive them from one ensemble member only, which might coincidentally be in a phase of high or low polynya activity.

The resulting polynya areas derived from observational and model sea ice thickness are not directly comparable for two reasons: the historical CMIP6 run (1 January 1850–31 December 2014) and the SMOS thin ice thickness data (1 June 2010 onwards) have less than 5 years of overlap. The second reason is their different definitions. The CMIP6 equivalent thickness variable describes the sea ice volume, or more exactly the sea ice mass per grid cell area (Notz et al., 2016), while the SMOS thin ice thickness describes the physical sea ice thickness (Huntemann et al., 2014).

4 Polynya statistics in CMIP6

The aim of our study is to determine whether the representation of polynyas in CMIP6 models is accurate in terms of location (Sect. 4.1) and frequency (Sect. 4.2). To do so, we compare the observational products and the model output with the same methods introduced in Sect. 3. We finally investigate the effect of the modelled polynyas on modelled stratification in Sect. 4.3.

4.1 Spatial distribution of Southern Ocean polynyas

We here evaluate the modelled location and spatial extent of coastal and open-water polynyas from the sea ice concentration and thickness variables. All 27 CMIP6 models exhibit coastal polynyas. The largest coastal polynya, the Ross Sea polynya, is represented in all models (Figs. 3, A3, A4), although some models open the Ross Sea polynya later than September (e.g. BCCs, CAMS-CSM1). Smaller polynyas, usually located by bays, headlands or islands along the continent, are best reproduced by models with a high horizontal resolution (Figs. 3, A3, A4). The GFDL-CM4 model with the highest horizontal resolution (ca. 0.25^∘) features the majority of polynyas seen in observations (Figs. 5 and A5). It is remarkable that all the coastal polynyas found in the observational product find their counterpart in the model. Even though GFDL-CM4 overestimates the polynya areas, the shapes and occurrence probability, for example of the Ross Sea polynya or the polynyas along the East Antarctic coast, are similar.

Figure 3 Spatial distribution of polynyas in the CMIP6 models, computed from the sea ice thickness output for September. The red-to-yellow colours indicate the number of years where polynyas occurred within the 165 years of the historical model run for each grid cell. The blue colours show the average sea ice concentration (in %).

In agreement with observations, the models show small coastal polynyas during winter that grow in size and in time, merging with other polynyas in spring/summer (November–January; see Figs. 3, 5b or A5). The coastal polynyas are generally located towards the west of geographical features along the continent outlines (Figs. 3, A3 and A4) where the polar easterlies prevail (Barber and Massom, 2007), except for polynyas at the southern end of the western Antarctic Peninsula, which is under the influence of the westerlies and the ACC and thus produces polynyas that are located towards the east of the land. Our analysis includes eight model families (e.g. ACCESS, BCC, CESM2…) that participate in CMIP6 with different model versions (Table 1). Seven of these families show polynyas at very similar locations for all members inside one family (Figs. 3, A3 and 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 (Figs. A4, A3) shows similar locations, but we found on average almost twice the polynya area with the sea ice thickness threshold method (Table 2).

Half of the 27 models show OWPs in the Ross or the Weddell seas. These polynyas are most common close to Maud Rise in the Weddell Sea (Figs. 3, A3 A4). Only 10 models show recurring OWPs in the Weddell Sea: ACCESS-ESM1-5, BCC-(CM2-MR/ESM1), CAMS-CSM1-0, EC-Earth3(-Veg), GFDL-ESM4, HadGEM3-GC31-LL and MPI-ESM1-2-(LR/HR). Of the 23 (28) models that provided sea ice thickness (concentration) output, 6 (8) show, on average, larger OWP than coastal polynya areas, while the remaining 17 (20) models show more coastal polynyas (Fig. 6). Five (fifteen) models overestimate the total polynya area compared to sea ice thickness (concentration) observations. Nine models (ACCESS-CM2, CAMS-CSM1.0, CNRM-CM6, CNRM-CSM2, GFDL-ESM4, all MPIs, UKESM) underestimate the polynya area when derived from the sea ice thickness (12 cm threshold) but at the same time overestimate it when derived from sea ice concentration (30 % threshold). Assuming the observational sea ice concentration product is accurate, this indicates these nine models show too low sea ice concentrations within regions covered by thin sea ice. Overall, the multi-model-spread in the polynya area is large, regardless of whether it is derived from monthly sea ice thickness, monthly concentration or daily concentration ranging from almost no polynyas (SAM0-UNICON) up to 214 930, 64 470 or 80 220 km² respectively for EC-Earth3-Veg. The multi-model mean is almost twice as high as the observed area (21 880 km² compared to 10 640 km² for monthly concentration, 38 560 to 12 050 km² for daily concentration and 58 680 to 73 060 km² for monthly sea ice thickness). This is driven primarily by five models; MPI-ESM1-2-HR, ACCESS-CM2, CNRM-ESM2 and the EC-Earth3s all overestimate the polynya area by a factor of 5. Surprisingly, the model ACCESS-CM2 shows only 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 Sect. 5.3.

Figure 4 (a–i) Unrealistic “polynya areas” computed with a sea ice thickness threshold method from the ice floe thickness output data. (j–k) SAM0-UNICON sea ice floe thickness and sea ice concentration for one time step. The red-to-yellow colours indicate the number of years where polynyas were detected by our algorithm within the 165 years of the historical model run for each grid cell. The blue colours show the average sea ice concentration (in %). The red contour line in panels (j)–(k) marks the 1 % sea ice concentration boundary. Open ocean and ice shelves both appear in grey.

Figure 5 Frequency and location of polynya occurrences in December (colours) with the positions of coastal polynyas (blue dots, locations from Arrigo and Van Dijken, 2003) for the common time interval of (a) the observed satellite sea ice (OSI450) and (b) the GFDL-CM4 CMIP6 model with the highest spatial resolution (25×25 km).

Figure 6 Winter polynya area for the CMIP6 models and the observational products (SMOS thin ice and OSI450). The polynyas are computed from (a) monthly sea ice concentration (b) daily sea ice concentration and (c) monthly sea ice thickness. The solid colours mark the yearly mean polynya area, while the transparent colours show the yearly maximum polynya area (Eqs. 1 and 2). The lower panels show the difference in yearly mean polynya area if derived from monthly/daily sea ice concentration (d) and monthly sea ice concentration/thickness (e), for the subset of models that provided both outputs. All data sets were trimmed to the winter period according to Eqs. (1) and (2) respectively, and then the yearly mean from 1979 to and including 2014 for homogenisation in time was formed; the only exception is the observational sea ice thickness, which is not available before 2010 and is presented in a separate panel at its full length until November 2020 at the side of panels (c) and (e). For the mean values in panels (a), (b) and (c), the black bars mark the range in which 50 % of the results are located.

Figure 7 Mean total winter polynya area as determined from monthly sea ice thickness for the CMIP6 models and the observational product (SMOS thin ice). For the models, results are shown for the full 165 years of the historical scenario, the observational data set includes the years 2010–2020. Each black dot is an average yearly polynya area for one model. The violin shape indicates the distribution of the values, the most values are in the thickest spot of the violins. The white dot indicates the mean value averaged over all available years.

4.2 Polynya frequency and seasonality

In Sect. 4.1, we discussed the average polynya areas over 35 years. A closer look at individual yearly values from the full historical run reveals that the polynya areas vary from year to year (Fig. 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 and highest for CESM2-FV2, HadGEM3-GC31-LL, MPI-ESM-1-2-HAM and SAM0-UNICON, which have 1 or more years without any coastal polynyas). Models with large coastal polynyas usually have little or no OWPs, while models with large OWPs have few coastal polynyas (Fig. 7). This is most distinct for comparison between the ACCESS models but can also bee seen in comparison of for example MPI to CESM2 models or CAMS-CSM1.0 to MRI-ESM2. The models with a larger fraction of coastal polynyas are ACCESS-CM2, all CESM2s, CNRMs, EC-Earth3s, IPSL-CM6A-LR, MRI-ESM2.0 and UKESM1-0-LL; for the rest the OWPs are dominating. Heuzé (2021) found a negative bottom density bias in the deep ocean for the latter, as they form the Antarctic Bottom Water frequently via deep convection instead of shelf processes. We believe that the deep convection is the reason for the OWPs, because the models that have OWPs mention convective overturning in the Southern Ocean in their model documentation (Table A1), and the water stratification within the polynyas indicates convection (Sect. 4.3). Models without the negative bottom density bias show predominantly coastal polynyas (CESM2s, CanESM5, CNRM-CM6.1, CanESM5, IPSL-CM6A, SAM0-UNICON, UKESM1-0-LL). In contrast to the stratification, numerical choices of parameterisation of open-ocean convective processes were not significant for the presence or absence of deep convection in CMIP5 models (De Lavergne et al., 2014).

The distribution of the OWPs has one fundamental difference to the one for the coastal polynyas. All models except for the BCCs show years without any or with negligible OWP activity (<10 000 km²). The OWP activity we have found in Sect. 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 (Figs. 8, A6 and A7). This is interesting since the sea ice in the Weddell Sea melts almost completely every summer. Thus polynya formation must either be facilitated by an external forcing that persists over several years (Campbell et al., 2019, suggest positive SAM fluctuations) or the vertical mixing in OWPs must act as a preconditioning for following years (Martinson et al., 1981). We will discuss this further in Sect. 4.3. While ACCESS-CM2 and BCC-(ESM1/CSM2) show OWPs almost every year, MPI-ESM1.2-(LR/HR), EC-Earth3 and CAMS-CSM1 show episodes of high OWP activity separated by multiple decades of little to no OWPs, as is the case for the real-world Weddell polynya (5 out of 50 years of observation, bottom of Figs. 8, A6, and A7 and Carsey, 1980; Campbell et al., 2019).

For the majority of models, the coastal polynyas appear randomly distributed over the years (Figs. 8, A6 and 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 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 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.

Figure 8 Yearly mean polynya area between May and November in the historical CMIP6 model run for the Weddell Sea region computed from the sea ice thickness. Each bar represents one CMIP6 model and is divided horizontally into the coastal polynya area (upper bar) and OWP area (lower bar). The values next to the model name represent the number of years an OWP with an area of more than 10.0×10³ km² was detected.

Figure 9 Seasonality of (a, d) the sea ice area; (b, e) the coastal polynyas area; and (c, f) the OWP area, for each CMIP6 model (symbols) and the observation products (plain black lines) in the Weddell Sea sector (65^∘ W–30^∘ E). The areas are derived from CMIP6 sea ice thickness and observational SMOS data (a, b, c) and from CMIP6 sea ice concentrations and observational OSI450 data (d, e, f). Grey shading indicates a 2 standard deviation range from the satellite sea ice data. The data were averaged for the time period 1979–2015.

4.3 Vertical ocean stratification in OWPs

OWPs are strongly influenced by upwelled sensible heat, and they are a source of bottom water production, both in observations (Cheon and Gordon, 2019) and in CMIP models (Heuzé et al., 2013); hence we expect to see a clear difference in the water stratification inside and outside an OWP.

We concentrate on the three models that form most OWPs (MPI-ESM1.2-HR, ACCESS-ESM1.5, BCC-ESM1; see Table 2) and show their vertical profiles in Fig. 10. For comparison, we present the observed hydrographic data of a SOCCOM profiling float (Johnson et al., 2018), which was deployed in January 2015 and surfaced two times in the Maud Rise polynya in winter 2017 (Campbell et al., 2019). To provide regionally and seasonally comparable data sets for the models and the profiling float, we chose to extract vertical profiles during the month of September from within a rectangle around the profiling float trajectory (see Campbell et al., 2019) with the edge coordinates 61–66^∘ S, 0–6^∘ E in the Weddell Sea. This region includes the northern flank of Maud Rise, where we found OWPs to be most common (Figs. 3, A3, A4).

For the SOCCOM profiling float, we use the information provided in Campbell et al. (2019) to differentiate vertical profiles when the float surfaces within an open-water polynya from those sampled under the sea ice. For the models, we use our algorithm to differentiate and group the grid points by whether they are within an OWP or not. Based on this criterion, we extract and group the monthly vertical salinity and conservative temperature profiles and average spatially over the different vertical profiles to obtain one averaged under-sea-ice profile and one averaged OWP profile (if OWP present in domain) per time step. We plot these profiles then in either blue (under sea ice) or red colour (OWP) in Fig. 10. In comparison to the CMIP6 models, the data of the SOCCOM float are of short length (2015–2017) and limited depth (maximum 2000 m). We hence use additional under-sea-ice profiles from the EN4.2.1 analyses full-depth ocean temperature and salinity data set (Good et al., 2013), limiting it to the same spatial domain as for the CMIP6 models, and show the profiles from 1980–2015 in light grey colour in Fig. 10m–p. We provide the EN4.2.1 data as additional reference of the observed water stratification without major OWP events

Table 2 Mean yearly areas of coastal polynyas (pa_co), OWPs (pa_op), their sum (pa_tot) and the maximum sea ice area (sa_tot) in 10³ km² derived from monthly/daily sea ice concentration and monthly sea ice thickness for September for the complete Southern Ocean, averaged over 1979–2015. Observational sea ice thickness marked with an asterisk (^*) is available over a 10-year period only. The four models matching the observations closest are printed in bold, and italic models have a an r value less than 0.8 (Figs. A1 and A2) and thus a relatively weak consistency over the three used sea ice variables.

In the three CMIP6 models, we find major differences in the vertical stratification between the profiles under sea ice and inside OWPs (Fig. 10a–l). In polynyas, deep mixing results in a very weakly stratified water column down to 4000 m, and in particular the salinity difference between mixed layer and Circumpolar Deep Water is reduced by up to 80 % (not shown). This is mainly due to an increased salinity in the upper 1000 m. Close inspection reveals that there is also a reduction in salinity in the Circumpolar Deep Water. Especially for the models MPI-ESM1.2-HR and ACCESS-ESM1.5 we also see a shoaling of the temperature maximum from 824 to 645 and 1560 to 1150 m depth, respectively, and increased temperature from 0.1 to 0.2 ^∘C for both models in the upper 500 m of the averaged profile. In the observational profiles (Fig. 10m–p), the temperature maximum can be found at a depth of about 200 m. Compared to the float data, ACCESS-ESM1.5 and BCC-ESM1 underestimate the temperature maximum by 1.5 and 1 ^∘C respectively and show it to occur too deep: 1148 ± 132 and 1366 ± 282 m, respectively, under the sea ice and 1560 ± 43 and 1205 ± 389 m, respectively, in OWPs.

In agreement with Cheon and Gordon (2019), we also find that the temperature profiles under the sea ice are more uniform than the profiles within the polynya (Fig. 10). The coldest profiles within polynyas are in fact the least stratified ones. For the model MPI-ESM1.2-HR (Fig. 10a–d), for example, the under-sea-ice temperature varies about 0.5 ^∘C from the mean value but varies by more than 1 ^∘C inside OWPs. Besides, inside OWPs, some of the profiles have a relatively shallow heat maximum (∼300 m), which is about 0.1 ^∘C warmer than the heat maximum under sea ice (Fig. 10b). A possible explanation could be that wind-driven upwelling contributed to a shallowing of the mixed layer depth (MLD) (hence the comparatively shallow heat maximum), leading to increased entrainment of warm CDW into the mixed layer and the opening of an OWP. Upon the opening of the polynya, the water cools rapidly while deep convection occurs, which could explain why the coldest profiles are the least stratified ones. The monthly resolution of the data sets is too low for us to verify these hypotheses.

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). Deep convection reaches different depths in the CMIP models (Fig. 10a–l). MPI-ESM1.2-HR contains some remaining stratification in all profiles, ACCESS-ESM1.5 shows evidence for mixing all the way to the bottom in many profiles and BCC-ESM1 conserves some stratification close to the bottom. In comparison to Argo float measurements, which we found to have undisturbed stable stratification below 500 m even in the polynya profiles, the models show deeper-reaching convection with a portion of the profiles showing no substantial stratification to the bottom (5000 m). However, not all profiles from inside an OWP are homogenised by mixing. About half of the profiles have a shallower and, for the MPI model, warmer temperature maximum compared to the under-ice case. This results in a shallower temperature maximum on average (Fig. 10a, e, i). In comparison to Argo float data, the vertical mixing in the models seems overestimated, which was to be expected as we here selected the subsample of three models with the highest OWP activity. Moreover, the float data contains only two profiles from the very edge of the 2017 Maud Rise polynya event (Campbell et al., 2019), where the vertical mixing may be less prominent than at the polynya's centre. A discussion of the stratification changes from the 1970s Weddell Sea polynya can be found in Cheon and Gordon (2019), where similar mixing as in the CMIP6 models was found.

Re-arranging the profiles into temperature–salinity diagrams (Fig. 10d, 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 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²) even under sea ice, where they hardly exceed $\mathrm{0.5}\times {\mathrm{10}}^{-\mathrm{5}}$ s⁻² (Fig. 10c, g, k, o). The deep mixing 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).

Figure 10 Winter profiles occurring within a polynya (red) and under the sea ice (blue, i.e. when no polynya is present) for three representative models (first three rows) and in observations (last row). From left to right: potential temperature, salinity, squared Brunt–Väisälä frequency (N²) and T–S diagram. Note the difference in vertical scale between the first two and the third column. The contour lines in the last column are lines of constant potential density. Two vertical Argo float profiles from within the polynya are available; for comparison we plot 19 under-sea-ice profiles (blue) and monthly September profiles from the EN4.2.1 climatology (grey) in panels (m)–(p).

5 Discussion

In the previous sections, we looked at the representation of polynyas in CMIP6 models and their effect on the local water stratification. We found large variations among models, with the majority underestimating OWP area but a few largely overestimating it. The 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 Current can explain its OWP representation. Then we examine the possible reasons for the strong overestimation of OWPs in some models and give an outlook about how we expect these biases to change in the near future.

5.1 Polynyas in CM versus ESM versions

The ACCESS, BCC and GFDL families each have two different model configurations for CMIP6: a climate model (CM) version and an Earth system model (ESM) version. In addition to climate physics, ESMs include physical, chemical, and biological processes and interactions and can therefore lead to more realistic climate predictions. Dong et al. (2020) found that the ESM versions have less climate sensitivity compared to the CM versions. We find that the ESM versions show more OWP activity (Fig. 6). The decreased OWP activity we find for CMs in our CMIP6 data set with ongoing global warming is consistent with the results of De Lavergne et al. (2014), in which OWPs in the Weddell Sea eventually stop at the end of the extended CMIP5 climate change runs, as the CM models show a stronger warming response than the ESM versions Dong et al. (2020).

CM models can usually be run at a higher resolution, due to their reduced complexity in comparison to ESMs. GFDL-CM4 has a higher oceanic resolution than GFDL-ESM4 and BCC-CSM2 has a higher atmospheric resolution than BCC-ESM1 (Table 1). A high spatial resolution can improve the representation of katabatic winds and coastlines, which might be why the above-mentioned CMs have a larger fraction of coastal polynyas than their ESM counterparts (Fig. 6). The CM model versions do not only show larger polynya areas (Figs. 6 and 3), but they also show more frequent occurrences of OWPs (Fig. 8). The same trend can be seen between the MPI models, which are alike except for their difference in resolution. The strongest difference in the positions and type of polynyas is found between the ACCESS-CM2 and ACCESS-ESM1.5 models, which show primarily coastal or open-water polynyas respectively (Figs. 3, A3, A4). This difference is further discussed in Sect. 5.3. Our comparison between CM and ESM was only possible on the 3 out of 15 model families that provided both types of model and as such may not apply to the other families. In contrast to our results, De Lavergne et al. (2014) found that deep convection in some CMIP5 models ceases earlier than in ESM models, but likewise their analysis was limited to only three families of models.

5.2 Connections between wind forcing, OWPs and the ACC

Open-water polynya activity is commonly attributed to the 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., 2016). We combine our polynya characteristics from Sect. 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) present 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.49: that is, the models with the strongest Southern Hemisphere westerlies have the highest OWP activity, while the models with the weakest westerlies show no OWP activity (Fig. 11a). We found this correlation to be the same for polynya detection from sea ice concentration and sea ice thickness. The relatively weak correlation coefficients reflects that OWPs cannot be attributed to the wind forcing alone but are sensitive to ocean and atmosphere preconditioning as well as the models' properties as mentioned in Sect. 1.

Figure 11 Connection of yearly OWP area (as in Eq. 2) with CMIP6 Southern Ocean properties from (Beadling et al., 2020). Total zonally averaged Southern Ocean wind forcing (x) and averaged OWP area (y) (a), linear regression plot of the maximum zonally integrated wind stress curl for each CMIP6 model and the latitude of that maximum (b), wind stress curl minimum (x) and OWP area (y) (c), and OWP area vs. ACC transport (d).

All observed Maud Rise polynyas were preceded by a positive SAM (Campbell et al., 2019). During a positive phase, the westerly wind belt intensifies and contracts towards Antarctica. This process leads to an increased cyclonic wind stress curl and intensifies the Weddell Gyre (Cheon et al., 2014). For the CMIP6 models we found that the maximum strength of the zonally averaged wind stress curl correlates with its latitude (Pearson correlation coefficient of −0.76, computed again using values from Beadling et al., 2020; see Fig. 11b). That means that, unsurprisingly, the CMIP6 models have different average SAM index values. Consequently, we expected an increased OWP activity in the Weddell Sea for the models with the highest wind stress curl maximum. This would consolidate the correlation between SAM and polynya activity. In contrast to Campbell et al. (2019) we find no significant correlation between the wind stress curl maximum and the average OWP area. CMIP6 models temporarily or constantly in a phase with a high SAM index do not show significantly increased OWP activity.

As shown by Hirabara et al. (2012), the ACC transport is positively correlated with OWP activity in observations. Models with an overestimation of OWPs also have too much deep convection (Heuzé, 2021, and Fig. 10), which increases the meridional density gradient and thus enhances the ACC transport (Hirabara et al., 2012; Beadling et al., 2020). The CMIP6 models that feature an ACC in agreement with observational values are ACCESS-CM2, ACCESS-ESM1-5, BCC-ESM1, CanESM5, GFDL-ESM4, IPSL-CM6A-LR, MIROC6, MPI-ESM-1-2-HAM, MPI-ESM1-2-LR and UKESM1-0-LL (Beadling et al., 2020). Nine of these 10 models exhibit open-ocean convection (Heuzé, 2021; see the models marked with an asterisk in Table A1), and in seven we found OWPs (Fig. 11d) in the Southern Ocean. The models ACCESS-CM2, MIROC6 and UKESM1-0-LL did not have large OWPs but excessive deep convection in overly large coastal polynyas (Heuzé, 2021, and Fig. 3) in the open Southern Ocean due to a strong negative sea ice bias (Fig. A3 and Heuzé, 2021; Tatebe et al., 2019) or periodical deep convection in the southern Weddell Sea (Heuzé, 2021). In summary, we find an across-model relationship between the strength of the ACC and OWP activity; models with an overestimation of open-ocean deep convection (which commonly leads to OWPs) are more likely to simulate a realistic ACC transport. Since most models simulate an ACC weaker than the observations (Beadling et al., 2020), an explanation could be that the models suffer in general from insufficient deep water formation, which is weakening the meridional density gradients and vertical stratification, but in turn lead to OWPs in some models that partially compensate for the lack of dense bottom water with spurious deep convection (Heuzé, 2015). Another explanation for this bias could be a general warm bias of the Circumpolar Deep Water (found for CMIP5 models Sallée et al., 2013) that causes OWPs rather than under-sea-ice deep convection (Dufour et al., 2017). We have shown that overabundant OWPs are a common issue in models (Sect. 4.1 and 4.2); that is a problem because it affects the whole water column (Sect. 4.3). We have now found relationships with the models' wider representation of the 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.

5.3 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 knowledge in the community given that within model documentation there are many descriptions of approaches to reduce deep convection and the presence of OWPs in the Southern Ocean (listed and highlighted with an asterisk in Table A1). We first describe approaches to reduce OWPs and conclude with future perspectives. In the previous section, we discussed the correlation of the wind stress curl and ACC transport to polynya activity; we could confirm this in a multi model correlation for the ACC. However, Beadling et al. (2020) found that out of the 34 CMIP6 models they analysed, 29 underestimated the ACC (of which only 12 significantly), 30 showed their wind stress curl not sufficiently south (5 significantly) and 30 underestimated the wind stress curl (WSC) minimum (9 significantly). All these parameters are positively correlated with OWP activity in observations (Campbell et al., 2019). If future model generations are approaching more realistic values for ACC strength, wind stress curl and SAM, we expect that the problem of open-water convection with the formation of large OWP will become evident for even more models.

One of the development foci for the GFDL models was to reduce Southern Ocean polynyas (Held et al., 2019). 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 (Held et al., 2019). One reason for this partial success only may be that a realistic representation of the Antarctic Slope Current and related Antarctic Slope Front restricts the lateral spread of fresh water away from the shelf, as Lockwood et al. (2021) found for the GFDL-CM2.6 climate model with a comparatively high horizontal resolution

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 find that the CMIP6 CESM2 models show coastal polynyas, but never OWPs. The observed polynya area is underestimated by at least 30 % by all analysed CESM2 models.

We found that both ACCESS models show more than double the polynya area of the observational data (Figs. 6 and 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 in the Weddell Sea form too often due to deep convection. An iceberg discharge scheme was hence introduced into the model ACCESS-CM2, which freshens and cools the upper ocean (Siobhan O'Farrell, personal communication, June 2020). Fresh water is transported further out into the open ocean instead of entering the ocean directly at the coast. This change effectively 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 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 examples of modellers' approaches to reduce deep convection and OWPs in the Southern Ocean. In the model documentations (listed in Table A1), we found more examples with the same aim. On the other hand, several models have found a way to prevent 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 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 Fig. 6). The Modular Ocean Model (MOM) introduced isopycnal coordinates in its latest version MOM6. Since several models are based on MOM, e.g. the ACCESS and the BCC models, we expect an improvement in deep water formation, Southern Ocean deep water properties and OWPs when these models successively upgrade to MOM6. Even though isopycnal coordinates were demonstrated to significantly reduce model drift and ocean heat uptake (Adcroft et al., 2019), they also have poor vertical resolution in weakly stratified regions such as the Weddell Sea and cannot avoid spurious deep convection in all cases. Only in combination with adequate implementation of shelf processes and dense water overflows, optimal results could be obtained (Adcroft et al., 2019).

6 Conclusions

In this paper, we evaluated the representation of Southern Ocean open-water and coastal polynyas in CMIP6 climate models and their effects on the modelled Weddell Sea. We found the following:

• all 27 analysed models have coastal polynyas around the Antarctic continent, while OWPs are present in only half of the models;

• CMIP6 models show OWPs most commonly either in the Weddell or the Ross seas;

• the position of polynya formation is very similar for models in one family and likely determined by the model properties;

• in comparison to observations, nine models underestimate polynya areas based on thickness threshold but overestimate them if based on the concentration threshold method;

• coastal polynya areas in CMIP6 have a large annual variability of at least a factor of 2.5;

• with total polynya areas from 6.5×10³ up to 215×10³ km², CMIP6 models show a large model spread.

Based on the results from the sea ice thickness, the models ACCESS-ESM1.5, BCC-ESM1 and MRI-ESM2.0 show the best agreement in total polynya area with less than 15 % deviation from observations, but ACCESS-ESM1.5 and BCC-ESM1 show a strong bias towards OWPs and MRI-ESM2.0 has only coastal polynyas. The inconsistent representation of polynyas can be problematic, since polynyas account for large amounts of sea ice and deep water production. We found larger polynyas with the sea ice thickness threshold method in 19 of 23 models and the observational data, and 4 models showed larger total polynya areas with the sea ice concentration threshold method. Using monthly sea ice concentration data instead of daily data leads to an underestimation of the polynya area in all models and the observations.

We have shown how OWPs significantly reduce the local stratification. We found indications of 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 Sect. 5.3 we discussed how the problem of spurious open-ocean deep convection in many CMIP6 models may be related to the OWPs biases. Half of the CMIP6 models show an over-representation of OWPs (Table 2), which is a known problem in the field (e.g. Held et al., 2019; Gutjahr et al., 2019; Sellar et al., 2019), likely caused by unrealistic deep convection events in the Southern Ocean (see Table A1, for references). We have described some of the models' strategies to prevent open-ocean polynyas and their caveats. CMIP6 models with realistic representation of the ACC transport and Southern Ocean wind fields commonly show unrealistically increased OWP formation in the Weddell Sea (Sect. 5.2). At the same time, future increases of the horizontal model resolution, improved vertical discretisation schemes, overflow parameterisations and freshwater influxes will likely help to achieve more realistic representations of polynyas in the Southern Ocean and in sea-ice-covered polar regions in general.

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 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 1 order of magnitude, we normalised 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 1 year; red are the coastal polynyas, blue the 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 1 order of magnitude, we normalised 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 A3 Spatial distribution of polynyas in the CMIP6 models, computed from the monthly sea ice concentration output for September. The red-to-yellow colours indicate the number of years where polynyas occurred within the 165 years of the historical model run for each grid cell. The blue colours show the average sea ice concentration (in %).

Figure A4 Spatial distribution of polynyas in the CMIP6 models, computed from the daily sea ice concentration output for September. The red-to-yellow colours indicate the number of years where polynyas occurred within the 165 years of the historical model run for each grid cell. The blue colours show the average sea ice concentration (in %).

Figure A5 The seasonal cycle of sea ice concentration (blue colours) and polynya probability (red–yellow) for the GFDL-CM4 model.

Figure A6 Yearly mean polynya area between May and November in the historical CMIP6 model run for the Weddell Sea region, computed from monthly sea ice concentrations. Each bar represents one CMIP6 model and is divided horizontally into coastal polynya area (upper bar) and OWP area (lower bar). The values next to the model name represent the number of years an open-water polynya with an area of more than 10.0×10³ km² was detected.

Figure A7 Yearly mean polynya area between May and November in the historical CMIP6 model run for the Weddell Sea region, computed from daily sea ice concentrations. Each bar represents one CMIP6 model and is divided horizontally into coastal polynya area (upper bar) and OWP area (lower bar). The values next to the model name represent the number of years an open-water polynya with an area of more than 10.0×10³ km² was detected.

Figure A8 Comparison of the winter profiles of the CMIP6 models ACCESS-ESM1.5 (blue) and ACCESS-CM2 (grey). From left to right: temperature, salinity, Brunt–Väisälä frequency N² and T–S diagram.

Bi et al. (2013)Ziehn et al. (2017)Bi et al. (2013)Wu et al. (2019)Wu et al. (2019)Rong et al. (2019)Danabasoglu et al. (2020)Danabasoglu et al. (2020)Danabasoglu et al. (2020)Danabasoglu et al. (2020)Voldoire et al. (2019)Séférian et al. (2019)Swart et al. (2019)Held et al. (2019)Dunne et al. (2019)Kelley et al. (2019)Kuhlbrodt et al. (2018)Boucher et al. (2020)Hajima et al. (2020)Tatebe et al. (2019)Mauritsen et al. (2019)Gutjahr et al. (2019)Müller et al. (2018)Mauritsen et al. (2019)Gutjahr et al. (2019)Yukimoto et al. (2019)Counillon et al. (2016)Park et al. (2019)Sellar et al. (2019)

Table A1 The 28 CMIP6 models used in this study (full name and references with model description). The references marked with an asterisk (^*) contain information about deep convection in the Southern Ocean. N/A indicates that no paper has been published yet for the CMIP6 configuration.