# Impacts of Vertical Convective Mixing Schemes and Freshwater Forcing on the 2016–2017 Maud Rise Polynya Openings in a Regional Ocean Simulation

Birte Gülk, Fabien Roquet, Alberto Naveira Garabato, Romain Bourdallé-Badie, Gurvan Madec, Hervé Giordani

<figure><img src="image_1.png" /></figure>

## To cite this version:

Birte Gülk, Fabien Roquet, Alberto Naveira Garabato, Romain Bourdallé-Badie, Gurvan Madec, et al. Impacts of Vertical Convective Mixing Schemes and Freshwater Forcing on the 2016–2017 Maud Rise Polynya Openings in a Regional Ocean Simulation. Journal of Advances in Modeling Earth Systems, 2024, 16 (5), 10.1029/2023ms004106. hal-04614188

HAL Id: hal-04614188

https://hal.science/hal-04614188

Submitted on 18 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Impacts of Vertical Convective Mixing Schemes and Freshwater Forcing on the 2016–2017 Maud Rise Polynya Openings in a Regional Ocean Simulation

Birte Gülk, Fabien Roquet, Alberto Naveira Garabato, Romain Bourdallé-Badie, Gurvan Madec, Hervé Giordani

► To cite this version:

Birte Gülk, Fabien Roquet, Alberto Naveira Garabato, Romain Bourdallé-Badie, Gurvan Madec, et al. Impacts of Vertical Convective Mixing Schemes and Freshwater Forcing on the 2016–2017 Maud Rise Polynya Openings in a Regional Ocean Simulation. Journal of Advances in Modeling Earth Systems, 2024, 16 (5), 10.1029/2023ms004106. hal-04614188

HAL Id: hal-04614188

https://hal.science/hal-04614188v1

Submitted on 18 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

10.1029/2023MS004106

**Key Points:**

* The Eddy-Diffusivity Mass-Flux (EDMF) parameterization is tested in a regional simulation of the ocean around Maud Rise

* Thermobaric effects on convective plumes are enabled by modifying the EDMF parameterization

* Simulations of Maud Rise polynya are highly sensitive to freshwater forcing and mixing schemes

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

B. Gülk,
birte.gulk@gu.se

**Citation:**

Gülk, B., Roquet, F., Naveira Garabato, A. C., Bourdallé-Badie, R., Madec, G., & Giordani, H. (2024). Impacts of vertical convective mixing schemes and freshwater forcing on the 2016–2017 Maud Rise polynya openings in a regional ocean simulation. *Journal of Advances in Modeling Earth Systems, 16*, e2023MS004106. https://doi.org/10.1029/2023MS004106

Received 6 NOV 2023

Accepted 4 APR 2024

© 2024 The Author(s). Journal of Advances in Modeling Earth Systems published by Wiley Periodicals LLC on behalf of American Geophysical Union. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_3.png" /></figure>

# Impacts of Vertical Convective Mixing Schemes and Freshwater Forcing on the 2016–2017 Maud Rise Polynya Openings in a Regional Ocean Simulation

Birte Gülk<fnref n="1" /> ID, Fabien Roquet<fnref n="1" /> ID, Alberto C. Naveira Garabato² ID, Romain Bourdallé-Badie³, Gurvan Madec<fnref n="4" />, and Hervé Giordani⁵ ID

**Abstract** The correct representation of the Maud Rise open-ocean polynya in the Weddell Sea remains a challenge for ocean models. Here we reproduce the most recent polynya openings in 2016–2017 using a regional configuration, and assess their dependencies on vertical convective mixing schemes and freshwater forcing, both separately and in combination. We test three vertical convective mixing schemes: the enhanced vertical diffusion (EVD), the Eddy-Diffusivity Mass-Flux (EDMF) parameterization, and a modified version of EDMF accounting for thermobaric effects. Using simulations for the period 2007–2017, we find that the modified EDMF reproduces the observed climatological evolution of the mixed layer depth better than the original EDMF and the EVD, but a polynya fails to open due to excessive freshwater forcing. We thus use the modified EDMF to perform sensitivity experiments with reduced precipitation during 2012–2017. The imposed freshwater forcing strongly affects the number of years with polynyas. The simulation with the best representation of the 2016–2017 polynyas is analyzed to evaluate the triggering mechanisms. The 2016 polynya was induced by the action of thermobaric instabilities on a weak ambient stratification. This opening preconditioned the water column for 2017, which produced a stronger polynya. By examining the impacts of the different convective mixing schemes, we show that the modified EDMF generates more realistic patterns of deep convection. Our results highlight the importance of surface freshwater forcing and thermobaricity in governing deep convection around Maud Rise, and the need to represent thermobaric instabilities to realistically model Maud Rise polynyas.

**Plain Language Summary** We investigate the impacts of representing numerical vertical mixing and surface freshwater forcing in a regional ocean model on polynyas (large openings in the pack ice) at Maud Rise, Southern Ocean. Maud Rise is prone to hosting polynyas, often associated with deep convection, which is a local vertical mixing process homogenizing the water column between surface and depths of several hundred meters. Numerical models often use simplistic strategies to represent this process, but improved parameterizations have recently become available. In this work, we test the impact of the representation of convective mixing in a particularly sensitive region. The last Maud Rise polynyas were observed in 2016 and 2017. Our regional simulation is capable of reproducing these polynyas, which has long been a challenge for ocean-sea ice models. We show that the 2016 polynya resulted from the action of a vertical instability at depth acting on weak ambient stratification. This event preconditioned the stronger 2017 polynya and deep convection. We conclude that representing convective plumes as a sub-grid scale process in models leads to a more realistic representation of open-ocean polynyas and associated convection events.

## 1. Introduction

The mechanisms of open-ocean polynya formation around Antarctica are poorly understood due to the complex nature of the interactions between ocean, sea ice and atmosphere that appear to govern polynya emergence. Open-ocean polynyas are offshore ice-free regions that occasionally emerge during winter. As a polynya opens, oceanic properties, such as heat, carbon dioxide and oxygen, are profoundly modified by exchanges with the atmosphere. The Weddell Sea is the Southern Ocean region with the strongest propensity to open-ocean polynya occurrence. There, Weddell Sea Polynyas (WSPs), large sea-ice openings with an area of $2–3 \cdot 10^5 \text{ km}^2$ were observed in 1974–1976 (Carsey, 1980). Over Maud Rise, smaller and short-lived polynyas, Maud Rise Polynyas (MRPs;

<footnote n="1">¹Department of Marine Sciences, University of Gothenburg, Gothenburg, Sweden, ²Ocean and Earth Science, University of Southampton, National Oceanography Centre, Southampton, UK, ³Mercator Ocean International, Toulouse, France,</footnote>

<footnote n="4">⁴LOCEAN-IPSL, Laboratoire d'Océanographie et du Climat: Expérimentation et Approches Numériques, Sorbonne Université/CNRS/IRD/MNHN, Paris, France, ⁵CNRM/CNRS, UMR-3589, Toulouse, France</footnote>
RESEARCH ARTICLE

10.1029/2023MS004106

**Key Points:**

* The Eddy-Diffusivity Mass-Flux (EDMF) parameterization is tested in a regional simulation of the ocean around Maud Rise

* Thermobaric effects on convective plumes are enabled by modifying the EDMF parameterization

* Simulations of Maud Rise polynya are highly sensitive to freshwater forcing and mixing schemes

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

B. Gülk,
birte.gulk@gu.se

**Citation:**

Gülk, B., Roquet, F., Naveira Garabato, A. C., Bourdallé-Badie, R., Madec, G., & Giordani, H. (2024). Impacts of vertical convective mixing schemes and freshwater forcing on the 2016–2017 Maud Rise polynya openings in a regional ocean simulation. *Journal of Advances in Modeling Earth Systems, 16*, e2023MS004106. https://doi.org/10.1029/2023MS004106

Received 6 NOV 2023

Accepted 4 APR 2024

© 2024 The Author(s). Journal of Advances in Modeling Earth Systems published by Wiley Periodicals LLC on behalf of American Geophysical Union. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_2.png" /></figure>

# Impacts of Vertical Convective Mixing Schemes and Freshwater Forcing on the 2016–2017 Maud Rise Polynya Openings in a Regional Ocean Simulation

Birte Gülk<fnref n="1" /> ID, Fabien Roquet<fnref n="1" /> ID, Alberto C. Naveira Garabato² ID, Romain Bourdallé-Badie³, Gurvan Madec<fnref n="4" />, and Hervé Giordani⁵ ID

**Abstract** The correct representation of the Maud Rise open-ocean polynya in the Weddell Sea remains a challenge for ocean models. Here we reproduce the most recent polynya openings in 2016–2017 using a regional configuration, and assess their dependencies on vertical convective mixing schemes and freshwater forcing, both separately and in combination. We test three vertical convective mixing schemes: the enhanced vertical diffusion (EVD), the Eddy-Diffusivity Mass-Flux (EDMF) parameterization, and a modified version of EDMF accounting for thermobaric effects. Using simulations for the period 2007–2017, we find that the modified EDMF reproduces the observed climatological evolution of the mixed layer depth better than the original EDMF and the EVD, but a polynya fails to open due to excessive freshwater forcing. We thus use the modified EDMF to perform sensitivity experiments with reduced precipitation during 2012–2017. The imposed freshwater forcing strongly affects the number of years with polynyas. The simulation with the best representation of the 2016–2017 polynyas is analyzed to evaluate the triggering mechanisms. The 2016 polynya was induced by the action of thermobaric instabilities on a weak ambient stratification. This opening preconditioned the water column for 2017, which produced a stronger polynya. By examining the impacts of the different convective mixing schemes, we show that the modified EDMF generates more realistic patterns of deep convection. Our results highlight the importance of surface freshwater forcing and thermobaricity in governing deep convection around Maud Rise, and the need to represent thermobaric instabilities to realistically model Maud Rise polynyas.

**Plain Language Summary** We investigate the impacts of representing numerical vertical mixing and surface freshwater forcing in a regional ocean model on polynyas (large openings in the pack ice) at Maud Rise, Southern Ocean. Maud Rise is prone to hosting polynyas, often associated with deep convection, which is a local vertical mixing process homogenizing the water column between surface and depths of several hundred meters. Numerical models often use simplistic strategies to represent this process, but improved parameterizations have recently become available. In this work, we test the impact of the representation of convective mixing in a particularly sensitive region. The last Maud Rise polynyas were observed in 2016 and 2017. Our regional simulation is capable of reproducing these polynyas, which has long been a challenge for ocean-sea ice models. We show that the 2016 polynya resulted from the action of a vertical instability at depth acting on weak ambient stratification. This event preconditioned the stronger 2017 polynya and deep convection. We conclude that representing convective plumes as a sub-grid scale process in models leads to a more realistic representation of open-ocean polynyas and associated convection events.

## 1. Introduction

The mechanisms of open-ocean polynya formation around Antarctica are poorly understood due to the complex nature of the interactions between ocean, sea ice and atmosphere that appear to govern polynya emergence. Open-ocean polynyas are offshore ice-free regions that occasionally emerge during winter. As a polynya opens, oceanic properties, such as heat, carbon dioxide and oxygen, are profoundly modified by exchanges with the atmosphere. The Weddell Sea is the Southern Ocean region with the strongest propensity to open-ocean polynya occurrence. There, Weddell Sea Polynyas (WSPs), large sea-ice openings with an area of $2–3 \cdot 10^5 \text{ km}^2$ were observed in 1974–1976 (Carsey, 1980). Over Maud Rise, smaller and short-lived polynyas, Maud Rise Polynyas (MRPs;

<footnote n="1">¹Department of Marine Sciences, University of Gothenburg, Gothenburg, Sweden, ²Ocean and Earth Science, University of Southampton, National Oceanography Centre, Southampton, UK, ³Mercator Ocean International, Toulouse, France,</footnote>

<footnote n="4">⁴LOCEAN-IPSL, Laboratoire d'Océanographie et du Climat: Expérimentation et Approches Numériques, Sorbonne Université/CNRS/IRD/MNHN, Paris, France, ⁵CNRM/CNRS, UMR-3589, Toulouse, France</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Topography of the Maud Rise region, with the schematized flow of the Halo (pink arrow) and the retroflection current (blue arrow). The green box indicates the area used to compute climatological variables. The blue (orange) line encircles the Taylor Cap (Flanks). The Taylor Cap is defined by the 2,500 m isobath, the Flanks are defined by the expansion of the 3,500 m isobath by ~55 km and limited on the northern flank by the 5,000 m isobath. The inlay shows the Weddell Sea and the zoomed-in area, which represents the model domain of Maud12 (red box).</figcaption></figure>

Comiso & Gordon, 1987), were last documented in 2016 and 2017. Maud Rise is a seamount rising from 5,000 m to 1,800 m at 2.5°E, 65°S, and is thought to localize the formation of WSPs and MRPs. The Rise stands in the way of the eastern limb of the Weddell Gyre, which conveys Warm Deep Water (WDW) westward into the Weddell Sea and partly feeds the Antarctic Slope Current (Leach et al., 2011). WDW originates from Circumpolar Deep Water (CDW) in the Antarctic Circumpolar Current. The Gyre's impingement on Maud Rise leads to the formation of a Taylor Cap atop the Rise, with a strong westward jet at the Rise's northern flank (the warm-water Halo) and a weaker jet at the southern flank (Gordon & Huber, 1990; Holland, 2001; Schröder & Fahrbach, 1999). An eastward retroflection current is also found at the Rise's northern flank, between the Taylor Cap and the Halo (Cisewski et al., 2011; Gülk et al., 2023). The warm-water Halo is reflected in the sea-ice properties, and leads to a thinning of the sea-ice concentration (SIC) above, while the Taylor Cap exhibits high SIC (Lindsay et al., 2004). See Figure 1 for a summary of the regional ocean dynamics and properties.

The opening of an open-ocean polynya is regularly accompanied by deep convection (Gordon, 1982), but deep convection can also occur without a polynya opening (Akitomo et al., 1995; Dufour et al., 2017). Deep convection homogenizes the properties between upper- and deep-ocean layers, and allows the deep ocean to exchange properties with the atmosphere. In the Weddell Sea, WDW/CDW has high carbon concentration due to its long-term isolation. When a convective event unfolds, carbon is released from the upwelled CDW. Depending on the air-sea CO₂ gradient and the event's duration, the CO₂ can partly be outgassed to the atmosphere or stored in the deep ocean (Bernardello et al., 2014). A deep convection event consists of several phases: the preconditioning phase, the deep convection phase, and the lateral exchange phase (Marshall & Schott, 1999). In the preconditioning phase, upper-ocean stratification is weakened, as a result of processes, for example, brine rejection or salt advection. The deep convection phase can be triggered by processes either at depth or at the surface, where such surface processes can be the same as those in the previous phase and increase upper-ocean density. At depth, thermobaric instabilities may destabilize the water column, inducing mixing of upper- and deep-ocean layers from below (Akitomo, 1999a). The deep convection phase results in the formation of dense convective plumes, which penetrate through the water column until they reach a denser layer. The lateral exchange phase can ensue simultaneously to the deep convective phase, and entails exchanges of properties between the convective plumes and the surrounding ocean (via e.g., de- or entrainment). Deep convection events are classified into two types (Akitomo, 1999a). Type I occurs in a weakly stratified ocean or within the mixed layer, and weakens the mixed
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Topography of the Maud Rise region, with the schematized flow of the Halo (pink arrow) and the retroflection current (blue arrow). The green box indicates the area used to compute climatological variables. The blue (orange) line encircles the Taylor Cap (Flanks). The Taylor Cap is defined by the 2,500 m isobath, the Flanks are defined by the expansion of the 3,500 m isobath by ~55 km and limited on the northern flank by the 5,000 m isobath. The inlay shows the Weddell Sea and the zoomed-in area, which represents the model domain of Maud12 (red box).</figcaption></figure>

Comiso & Gordon, 1987), were last documented in 2016 and 2017. Maud Rise is a seamount rising from 5,000 m to 1,800 m at 2.5°E, 65°S, and is thought to localize the formation of WSPs and MRPs. The Rise stands in the way of the eastern limb of the Weddell Gyre, which conveys Warm Deep Water (WDW) westward into the Weddell Sea and partly feeds the Antarctic Slope Current (Leach et al., 2011). WDW originates from Circumpolar Deep Water (CDW) in the Antarctic Circumpolar Current. The Gyre's impingement on Maud Rise leads to the formation of a Taylor Cap atop the Rise, with a strong westward jet at the Rise's northern flank (the warm-water Halo) and a weaker jet at the southern flank (Gordon & Huber, 1990; Holland, 2001; Schröder & Fahrbach, 1999). An eastward retroflection current is also found at the Rise's northern flank, between the Taylor Cap and the Halo (Cisewski et al., 2011; Gülk et al., 2023). The warm-water Halo is reflected in the sea-ice properties, and leads to a thinning of the sea-ice concentration (SIC) above, while the Taylor Cap exhibits high SIC (Lindsay et al., 2004). See Figure 1 for a summary of the regional ocean dynamics and properties.

The opening of an open-ocean polynya is regularly accompanied by deep convection (Gordon, 1982), but deep convection can also occur without a polynya opening (Akitomo et al., 1995; Dufour et al., 2017). Deep convection homogenizes the properties between upper- and deep-ocean layers, and allows the deep ocean to exchange properties with the atmosphere. In the Weddell Sea, WDW/CDW has high carbon concentration due to its long-term isolation. When a convective event unfolds, carbon is released from the upwelled CDW. Depending on the air-sea CO₂ gradient and the event's duration, the CO₂ can partly be outgassed to the atmosphere or stored in the deep ocean (Bernardello et al., 2014). A deep convection event consists of several phases: the preconditioning phase, the deep convection phase, and the lateral exchange phase (Marshall & Schott, 1999). In the preconditioning phase, upper-ocean stratification is weakened, as a result of processes, for example, brine rejection or salt advection. The deep convection phase can be triggered by processes either at depth or at the surface, where such surface processes can be the same as those in the previous phase and increase upper-ocean density. At depth, thermobaric instabilities may destabilize the water column, inducing mixing of upper- and deep-ocean layers from below (Akitomo, 1999a). The deep convection phase results in the formation of dense convective plumes, which penetrate through the water column until they reach a denser layer. The lateral exchange phase can ensue simultaneously to the deep convective phase, and entails exchanges of properties between the convective plumes and the surrounding ocean (via e.g., de- or entrainment). Deep convection events are classified into two types (Akitomo, 1999a). Type I occurs in a weakly stratified ocean or within the mixed layer, and weakens the mixed
layer base gradually. Type II occurs in a layered system, in which a cold and fresh layer lies above a warmer and more saline layer. Here, thermobaricity can cause an abrupt overturning of the water column (Akitomo, 1999a). Further, Type I may occur annually, while Type II takes place more irregularly, depending on ambient stratification conditions. A common region for Type I convection is the Greenland basin, while Type II can be found in the Weddell Sea (Akitomo, 1999b).

Several theories have been put forward regarding the triggering of WSPs and MRPs. Atmospheric observations suggest that MRP formation in 2016–2017 could have been initiated by polar cyclones (Francis et al., 2019) or by atmospheric rivers (Francis et al., 2020). In contrast, ocean observations point to an increase in the mixed layer salinity before the opening of a polynya by severe storms (Campbell et al., 2019). Observations have also shown that the Maud Rise region is susceptible to thermobaricity, which can lead to abrupt deep convection before a polynya opening (Akitomo, 2006; McPhee, 2003). Satellite observations revealed a thinning of sea ice 4 months ahead of the polynya opening caused by the atmospheric and oceanic circulation (Zhou et al., 2022). Other modeling and observational studies have emphasized the potential polynya-triggering role of eddy generation at the flanks of Maud Rise, resulting from the Weddell Gyre’s impingement on the Rise (Cheon & Gordon, 2019; Holland, 2001). As well as an increased subsurface temperature maximum atop Maud Rise, which preceded the latest openings (Gülk et al., 2023), climate-scale ocean models have suggested that the advection of a surface salinity anomaly from Astrid Ridge to Maud Rise (Kurtakoti et al., 2018) and the subsurface build-up of heat (Dufour et al., 2017; Rheinlænder et al., 2021) may have promoted polynya formation.

In summary, many processes may play a significant role in open-ocean polynya formation and deep convection around Maud Rise, but the overall mechanistic picture of these phenomena remains complex and incomplete. Such incomplete understanding stems from the scarcity of long-term observations, as well as the difficulties with reproducing reasonably realistic polynyas and deep convection events in numerical models. Earth System Models may reproduce polynyas (Kurtakoti et al., 2018; Rheinlænder et al., 2021), and be useful in studying the atmosphere-ocean-ice processes at play. A caveat of these models, however, is that the scales of polynyas (e.g., size, location and time) are often overestimated, such that polynya properties cannot easily translate to real-world scenarios. The Earth System Models included in CMIP6, tend to generate either no open-ocean polynya, or too frequent and too large polynyas (Mohrmann et al., 2021). Some ocean general circulation models (OGCMs) forced by atmospheric reanalyzes have been able to reproduce the 1970s WSP (Cheon et al., 2015), but not the smaller and short-lived MRPs. OGCMs with data assimilation, such as the Southern Ocean State Estimate (Mazloff et al., 2010) and GLORYS12 (Lellouche et al., 2021) exhibit a polynya, but the role of assimilation can be unclear. A major issue with polynya formation in OGCMs is the potentially acute sensitivity of their representation of polynya formation and deep convection events to a range of forcings, model parameters and resolved processes. For example, some studies have shown that relatively modest changes in the freshwater balance (initial sea-ice cover and surface forcing) in atmospherically-forced ocean-sea ice models of the Southern Ocean can substantially increase or reduce deep convection activity (Kjellsson et al., 2015; Stoessel et al., 2015). Further, vertical mixing parameters can hamper or amplify the occurrence of deep convection in polynya regions (Heuzé et al., 2015; Kjellsson et al., 2015). Important processes such as thermobaricity, which are thought to play a significative role in triggering deep convection events (Akitomo, 2006; McPhee, 2003), are not correctly represented in ocean models.

Understanding how open-ocean polynyas form around Maud Rise is a necessary step in advancing our ability to predict future openings of WSPs, which grow out of MRPs, and their consequences for Southern Ocean deep-water formation (Kurtakoti et al., 2021). Here, we seek to address this knowledge gap by investigating the polynya formation impacts of different vertical convective mixing schemes and the surface freshwater forcing in a regional ocean simulation encompassing the 2016–2017 MRP period and the preceding period from 2007 onward. For this purpose, we use a regional configuration of the Nucleus for European Modeling of the Ocean (NEMO) version 4, simulating the ocean and sea-ice conditions and perform a number of sensitivity runs. Our goal is to identify the conditions needed to reproduce MRPs and deep convection events. The model setup, the vertical convective mixing schemes and the sensitivity runs are presented in detail in Section 2.1. The observational data sets and methods used are described in Sections 2.2 and 2.3, and our results in Section 3. In Section 4, our main findings are discussed, and conclusions are drawn.
layer base gradually. Type II occurs in a layered system, in which a cold and fresh layer lies above a warmer and more saline layer. Here, thermobaricity can cause an abrupt overturning of the water column (Akitomo, 1999a). Further, Type I may occur annually, while Type II takes place more irregularly, depending on ambient stratification conditions. A common region for Type I convection is the Greenland basin, while Type II can be found in the Weddell Sea (Akitomo, 1999b).

Several theories have been put forward regarding the triggering of WSPs and MRPs. Atmospheric observations suggest that MRP formation in 2016–2017 could have been initiated by polar cyclones (Francis et al., 2019) or by atmospheric rivers (Francis et al., 2020). In contrast, ocean observations point to an increase in the mixed layer salinity before the opening of a polynya by severe storms (Campbell et al., 2019). Observations have also shown that the Maud Rise region is susceptible to thermobaricity, which can lead to abrupt deep convection before a polynya opening (Akitomo, 2006; McPhee, 2003). Satellite observations revealed a thinning of sea ice 4 months ahead of the polynya opening caused by the atmospheric and oceanic circulation (Zhou et al., 2022). Other modeling and observational studies have emphasized the potential polynya-triggering role of eddy generation at the flanks of Maud Rise, resulting from the Weddell Gyre’s impingement on the Rise (Cheon & Gordon, 2019; Holland, 2001). As well as an increased subsurface temperature maximum atop Maud Rise, which preceded the latest openings (Gülk et al., 2023), climate-scale ocean models have suggested that the advection of a surface salinity anomaly from Astrid Ridge to Maud Rise (Kurtakoti et al., 2018) and the subsurface build-up of heat (Dufour et al., 2017; Rheinlænder et al., 2021) may have promoted polynya formation.

In summary, many processes may play a significant role in open-ocean polynya formation and deep convection around Maud Rise, but the overall mechanistic picture of these phenomena remains complex and incomplete. Such incomplete understanding stems from the scarcity of long-term observations, as well as the difficulties with reproducing reasonably realistic polynyas and deep convection events in numerical models. Earth System Models may reproduce polynyas (Kurtakoti et al., 2018; Rheinlænder et al., 2021), and be useful in studying the atmosphere-ocean-ice processes at play. A caveat of these models, however, is that the scales of polynyas (e.g., size, location and time) are often overestimated, such that polynya properties cannot easily translate to real-world scenarios. The Earth System Models included in CMIP6, tend to generate either no open-ocean polynya, or too frequent and too large polynyas (Mohrmann et al., 2021). Some ocean general circulation models (OGCMs) forced by atmospheric reanalyzes have been able to reproduce the 1970s WSP (Cheon et al., 2015), but not the smaller and short-lived MRPs. OGCMs with data assimilation, such as the Southern Ocean State Estimate (Mazloff et al., 2010) and GLORYS12 (Lellouche et al., 2021) exhibit a polynya, but the role of assimilation can be unclear. A major issue with polynya formation in OGCMs is the potentially acute sensitivity of their representation of polynya formation and deep convection events to a range of forcings, model parameters and resolved processes. For example, some studies have shown that relatively modest changes in the freshwater balance (initial sea-ice cover and surface forcing) in atmospherically-forced ocean-sea ice models of the Southern Ocean can substantially increase or reduce deep convection activity (Kjellsson et al., 2015; Stoessel et al., 2015). Further, vertical mixing parameters can hamper or amplify the occurrence of deep convection in polynya regions (Heuzé et al., 2015; Kjellsson et al., 2015). Important processes such as thermobaricity, which are thought to play a significative role in triggering deep convection events (Akitomo, 2006; McPhee, 2003), are not correctly represented in ocean models.

Understanding how open-ocean polynyas form around Maud Rise is a necessary step in advancing our ability to predict future openings of WSPs, which grow out of MRPs, and their consequences for Southern Ocean deep-water formation (Kurtakoti et al., 2021). Here, we seek to address this knowledge gap by investigating the polynya formation impacts of different vertical convective mixing schemes and the surface freshwater forcing in a regional ocean simulation encompassing the 2016–2017 MRP period and the preceding period from 2007 onward. For this purpose, we use a regional configuration of the Nucleus for European Modeling of the Ocean (NEMO) version 4, simulating the ocean and sea-ice conditions and perform a number of sensitivity runs. Our goal is to identify the conditions needed to reproduce MRPs and deep convection events. The model setup, the vertical convective mixing schemes and the sensitivity runs are presented in detail in Section 2.1. The observational data sets and methods used are described in Sections 2.2 and 2.3, and our results in Section 3. In Section 4, our main findings are discussed, and conclusions are drawn.
## 2. Data and Methods

### 2.1. Model Set-Up and Experiments

We use a 1/12° horizontal resolution (equivalent to 3.9 km at 65°S) regional NEMO configuration, Maud12, with ocean (version 4.2.0) and sea ice components (SI3) (Madec & NEMO System Team, 2019; NEMO Sea Ice Working Group, 2018). The domain covers the region 5°W to 19.5°E and 61.5°S to 70.5°S (Figure 1 red box), and provides daily output on 50 vertical levels. The vertical levels increase in spacing with depth. At the surface a grid resolution of 1 m is used, and this increases to 500 m at the bottom. The atmospheric forcing is derived from JRA55-do version 1.5 (Tsujino et al., 2018). The open boundary conditions for the ocean and ice properties are derived from GLORYS12 (Lellouche et al., 2021). More details about the configuration can be found in Gülk et al. (2023). Adopting NEMO 4.2.0 allows the use of the Eddy-Diffusivity Mass-Flux (EDMF) parameterization to represent convection in the ocean. This parameterization has so far been tested and validated in 1-D scenarios within the Gulf of Lion (Mediterranean Sea) and at the PAPA station (North Pacific Ocean) (Giordani et al., 2020).

With Maud12, we first assess the impacts of the commonly used enhanced vertical diffusion (EVD) scheme, the EDMF scheme, and the modifications applied to a version of the EDMF scheme, which allow for thermobaric instabilities. The EDMF scheme and its modifications are described in Section 2.1.1. The schemes are tested in simulations spanning the period 2007–2017. These runs are used to investigate how the schemes affect the climatological evolution of the temperature, salinity and the mixed layer depth (MLD), and elicit caveats in the simulations. Second, we perform sensitivity runs with different surface freshwater forcings, using the best-performing convective mixing scheme from the first set of tests (Section 2.1.2). In these sensitivity runs, we focus on the occurrence of polynya openings and identify which freshwater forcing produces results in agreement with satellite observations. We further determine the processes governing the openings in 2016 and 2017. Using the optimal freshwater forcing we, again, test the different vertical convective mixing schemes, to assess the dependence of the key model results on the chosen scheme.

#### 2.1.1. Vertical Convective Mixing Schemes

The difference between the EVD and EDMF schemes resides in the way they reproduce convective plumes in numerical models. If the stratification is unstable ($N^2 < 0$), EVD assigns large eddy mixing coefficients (Madec & NEMO System Team, 2019) and thereby enhances the vertical mixing across two neighboring grid cells. In contrast, EDMF produces sub-grid scale convective plumes penetrating through the water column, by combining local turbulence induced by eddies and mass fluxes of a larger scale. The convective plumes are active in regions where the surface buoyancy flux is positive (i.e., out of the ocean). In stratified zones, the plume is accelerated/decelerated by potential density differences between the plume and surrounding waters. The modified version of EDMF contains several changes to the original scheme, and enables thermobaric effects to trigger convective plumes at depth by using in situ density instead of potential density to compute buoyancy differences for density variation.

In the EDMF scheme as presented by Giordani et al. (2020), the surface buoyancy flux was computed with fixed values for the thermal expansion ($\alpha$) and haline contraction ($\beta$) coefficients of $\alpha = 2 \cdot 10^{-4} \text{ K}^{-1}$ and $\beta = 7 \cdot 10^{-4} \text{ kg g}^{-1}$. The surface buoyancy flux is used to derive the fractional convective area of the plume, from which the plume is initiated. The plume velocity ($w_p$) is computed by integrating the equation,

$$ \left(\frac{1}{2} + a_2\right) \frac{\partial w_p^2}{\partial z} = a_1 F_b - a_2 g \frac{\rho_p}{\rho_h} w_p^2 \quad (1) $$

using the gravitational acceleration $g$, the plume buoyancy acceleration $F_b$, the surface-referenced potential density of the plume $\rho_p$, the hydrostatic pressure $p_h$ and the coefficients $a_1$ and $a_2$. Note that Giordani et al. (2020) used notation $\alpha$ for $a_2$, but we changed it to avoid confusion with the coefficient of thermal expansion. The value $a_1 = 1$ was diagnosed from Large Eddy Simulations in the atmosphere (Siebesma et al., 2003) and $a_2 = 0.2$ produced the best results in the EDMF 1-D test cases. $a_2 g \frac{\rho_p}{\rho_h}$ represents the resistance term associated with the non-hydrostatic pressure perturbation applied to downward parcels. As Equation 1 indicates the main driver of the
## 2. Data and Methods

### 2.1. Model Set-Up and Experiments

We use a 1/12° horizontal resolution (equivalent to 3.9 km at 65°S) regional NEMO configuration, Maud12, with ocean (version 4.2.0) and sea ice components (SI3) (Madec & NEMO System Team, 2019; NEMO Sea Ice Working Group, 2018). The domain covers the region 5°W to 19.5°E and 61.5°S to 70.5°S (Figure 1 red box), and provides daily output on 50 vertical levels. The vertical levels increase in spacing with depth. At the surface a grid resolution of 1 m is used, and this increases to 500 m at the bottom. The atmospheric forcing is derived from JRA55-do version 1.5 (Tsujino et al., 2018). The open boundary conditions for the ocean and ice properties are derived from GLORYS12 (Lellouche et al., 2021). More details about the configuration can be found in Gülk et al. (2023). Adopting NEMO 4.2.0 allows the use of the Eddy-Diffusivity Mass-Flux (EDMF) parameterization to represent convection in the ocean. This parameterization has so far been tested and validated in 1-D scenarios within the Gulf of Lion (Mediterranean Sea) and at the PAPA station (North Pacific Ocean) (Giordani et al., 2020).

With Maud12, we first assess the impacts of the commonly used enhanced vertical diffusion (EVD) scheme, the EDMF scheme, and the modifications applied to a version of the EDMF scheme, which allow for thermobaric instabilities. The EDMF scheme and its modifications are described in Section 2.1.1. The schemes are tested in simulations spanning the period 2007–2017. These runs are used to investigate how the schemes affect the climatological evolution of the temperature, salinity and the mixed layer depth (MLD), and elicit caveats in the simulations. Second, we perform sensitivity runs with different surface freshwater forcings, using the best-performing convective mixing scheme from the first set of tests (Section 2.1.2). In these sensitivity runs, we focus on the occurrence of polynya openings and identify which freshwater forcing produces results in agreement with satellite observations. We further determine the processes governing the openings in 2016 and 2017. Using the optimal freshwater forcing we, again, test the different vertical convective mixing schemes, to assess the dependence of the key model results on the chosen scheme.

#### 2.1.1. Vertical Convective Mixing Schemes

The difference between the EVD and EDMF schemes resides in the way they reproduce convective plumes in numerical models. If the stratification is unstable ($N^2 < 0$), EVD assigns large eddy mixing coefficients (Madec & NEMO System Team, 2019) and thereby enhances the vertical mixing across two neighboring grid cells. In contrast, EDMF produces sub-grid scale convective plumes penetrating through the water column, by combining local turbulence induced by eddies and mass fluxes of a larger scale. The convective plumes are active in regions where the surface buoyancy flux is positive (i.e., out of the ocean). In stratified zones, the plume is accelerated/decelerated by potential density differences between the plume and surrounding waters. The modified version of EDMF contains several changes to the original scheme, and enables thermobaric effects to trigger convective plumes at depth by using in situ density instead of potential density to compute buoyancy differences for density variation.

In the EDMF scheme as presented by Giordani et al. (2020), the surface buoyancy flux was computed with fixed values for the thermal expansion ($\alpha$) and haline contraction ($\beta$) coefficients of $\alpha = 2 \cdot 10^{-4} \text{ K}^{-1}$ and $\beta = 7 \cdot 10^{-4} \text{ kg g}^{-1}$. The surface buoyancy flux is used to derive the fractional convective area of the plume, from which the plume is initiated. The plume velocity ($w_p$) is computed by integrating the equation,

$$ \left(\frac{1}{2} + a_2\right) \frac{\partial w_p^2}{\partial z} = a_1 F_b - a_2 g \frac{\rho_p}{\rho_h} w_p^2 \quad (1) $$

using the gravitational acceleration $g$, the plume buoyancy acceleration $F_b$, the surface-referenced potential density of the plume $\rho_p$, the hydrostatic pressure $p_h$ and the coefficients $a_1$ and $a_2$. Note that Giordani et al. (2020) used notation $\alpha$ for $a_2$, but we changed it to avoid confusion with the coefficient of thermal expansion. The value $a_1 = 1$ was diagnosed from Large Eddy Simulations in the atmosphere (Siebesma et al., 2003) and $a_2 = 0.2$ produced the best results in the EDMF 1-D test cases. $a_2 g \frac{\rho_p}{\rho_h}$ represents the resistance term associated with the non-hydrostatic pressure perturbation applied to downward parcels. As Equation 1 indicates the main driver of the
plume's velocity is the buoyancy acceleration $F_b$, which is derived from the surface-referenced potential density difference between the plume ($\rho_p$) and its environment ($\rho_e$). The plume buoyancy acceleration is estimated by

$$F_b = g \frac{\rho_p - \rho_e}{\rho_e} \qquad (2)$$

When $F_b > 0$ ($F_b < 0$) the plume accelerates (decelerates). The buoyancy acceleration may lead to an overshoot below the maximum depth of static instability.

We have made several successive modifications on the EDMF scheme. In a first modified version, called EDMF$_{α}$, we used local values of $\alpha$ and $\beta$ as function of conservative temperature and absolute salinity for the surface buoyancy flux. This values are obtained using the TEOS-10 standard for the equation of state (Roquet et al., 2015). In the EDMF$_{ρ}$ version, we also modified the plume buoyancy acceleration (Equation 2) to use in situ density instead of potential density, again based on TEOS-10, enabling us to capture thermobaric effects. Finally, in the last modified version (EDMF$_{e}$), we added a frictional term $\epsilon_\psi w_p^2$ to the plume velocity equation (Equation 1) to include lateral entrainment of momentum (Gregory, 2001), and we used the hydrostatic approximation $g\rho_p/p_h \approx 1/h$ to compute the first frictional term, where $h$ is the vertical depth.

The plume velocity equation used in EDMF$_{e}$ becomes,

$$\left(\frac{1}{2} + a_2\right) \frac{\partial w_p^2}{\partial z} = a_1 F_b - \left(\frac{a_2}{h} + \epsilon_\psi\right) w_p^2 \qquad (3)$$

where the lateral entrainment rate $\epsilon_\psi$ is computed as follows,

$$\epsilon_\psi = \left| \frac{1}{w_p} \frac{\partial w_p}{\partial z} \right| + \epsilon_g \qquad (4)$$

The lateral entrainment rate depends on the small-scale turbulence at the plume-environment interface (first left hand side term) added to a background entrainment rate $\epsilon_g$. We used a constant value $\epsilon_g = 0.001$ throughout the water column.

### 2.1.2. Sensitivity Experiments

We generate two sets of sensitivity experiments: one with perturbed surface freshwater forcing using only one of the mixing schemes (EDMF$_{e}$), and one with each of the convective mixing schemes presented above using the same surface forcing. All sensitivity runs are branching off the EDMF$_{e}$ with unchanged forcing after 2011 (except for the surface freshwater flux in the respective experiments). This spin up simulation is chosen due to its best representation of oceanic properties amongst all runs conducted in the first set of tests (Section 3.1). An overview of all conducted runs is provided in Table 1.

In the sensitivity runs with perturbed freshwater forcing, the rain and snow precipitation is reduced in 2012–2017, whilst using EDMF$_{e}$. The precipitation is sequentially decreased to 0%, 10%, 20%, 30% and 40% of the original precipitation. A sequential decrease is necessary due to uncertainties in the freshwater import from the open-boundary sea-ice conditions and the precipitation in atmospheric reanalysis product. A comparison of the total precipitation between the used JRA55-do version 1.5 and ERA5 (Hersbach et al., 2020) shows a spatial variability (Figure S1 in Supporting Information S1). The sensitivity runs are validated by reference to the trend of the freshwater content in the domain (as explained in Section 2.3) as well as the polynya properties. Polynya properties are compared to SIC from satellite observations. The run with highest coherence with observations is then used to mechanically investigate polynya formation and deep convection events.

In the sensitivity runs involving different convective mixing schemes, we use the freshwater forcing from the run that exhibited the most realistic polynya properties, and replaced the convective mixing scheme EDMF$_{e}$ with the following ones: EVD, EDMF, EDMF$_{α}$, and EDMF$_{ρ}$. This test is expected to highlight the dependence of our results to the selected scheme. Additionally, simulations with 10% and 0% of the original precipitation using EVD
plume's velocity is the buoyancy acceleration $F_b$, which is derived from the surface-referenced potential density difference between the plume ($\rho_p$) and its environment ($\rho_e$). The plume buoyancy acceleration is estimated by

$$F_b = g \frac{\rho_p - \rho_e}{\rho_e} \qquad (2)$$

When $F_b > 0$ ($F_b < 0$) the plume accelerates (decelerates). The buoyancy acceleration may lead to an overshoot below the maximum depth of static instability.

We have made several successive modifications on the EDMF scheme. In a first modified version, called EDMF$_{α}$, we used local values of $\alpha$ and $\beta$ as function of conservative temperature and absolute salinity for the surface buoyancy flux. This values are obtained using the TEOS-10 standard for the equation of state (Roquet et al., 2015). In the EDMF$_{ρ}$ version, we also modified the plume buoyancy acceleration (Equation 2) to use in situ density instead of potential density, again based on TEOS-10, enabling us to capture thermobaric effects. Finally, in the last modified version (EDMF$_{e}$), we added a frictional term $\epsilon_\psi w_p^2$ to the plume velocity equation (Equation 1) to include lateral entrainment of momentum (Gregory, 2001), and we used the hydrostatic approximation $g\rho_p/p_h \approx 1/h$ to compute the first frictional term, where $h$ is the vertical depth.

The plume velocity equation used in EDMF$_{e}$ becomes,

$$\left(\frac{1}{2} + a_2\right) \frac{\partial w_p^2}{\partial z} = a_1 F_b - \left(\frac{a_2}{h} + \epsilon_\psi\right) w_p^2 \qquad (3)$$

where the lateral entrainment rate $\epsilon_\psi$ is computed as follows,

$$\epsilon_\psi = \left| \frac{1}{w_p} \frac{\partial w_p}{\partial z} \right| + \epsilon_g \qquad (4)$$

The lateral entrainment rate depends on the small-scale turbulence at the plume-environment interface (first left hand side term) added to a background entrainment rate $\epsilon_g$. We used a constant value $\epsilon_g = 0.001$ throughout the water column.

### 2.1.2. Sensitivity Experiments

We generate two sets of sensitivity experiments: one with perturbed surface freshwater forcing using only one of the mixing schemes (EDMF$_{e}$), and one with each of the convective mixing schemes presented above using the same surface forcing. All sensitivity runs are branching off the EDMF$_{e}$ with unchanged forcing after 2011 (except for the surface freshwater flux in the respective experiments). This spin up simulation is chosen due to its best representation of oceanic properties amongst all runs conducted in the first set of tests (Section 3.1). An overview of all conducted runs is provided in Table 1.

In the sensitivity runs with perturbed freshwater forcing, the rain and snow precipitation is reduced in 2012–2017, whilst using EDMF$_{e}$. The precipitation is sequentially decreased to 0%, 10%, 20%, 30% and 40% of the original precipitation. A sequential decrease is necessary due to uncertainties in the freshwater import from the open-boundary sea-ice conditions and the precipitation in atmospheric reanalysis product. A comparison of the total precipitation between the used JRA55-do version 1.5 and ERA5 (Hersbach et al., 2020) shows a spatial variability (Figure S1 in Supporting Information S1). The sensitivity runs are validated by reference to the trend of the freshwater content in the domain (as explained in Section 2.3) as well as the polynya properties. Polynya properties are compared to SIC from satellite observations. The run with highest coherence with observations is then used to mechanistically investigate polynya formation and deep convection events.

In the sensitivity runs involving different convective mixing schemes, we use the freshwater forcing from the run that exhibited the most realistic polynya properties, and replaced the convective mixing scheme EDMF$_{e}$ with the following ones: EVD, EDMF, EDMF$_{α}$, and EDMF$_{ρ}$. This test is expected to highlight the dependence of our results to the selected scheme. Additionally, simulations with 10% and 0% of the original precipitation using EVD
<table><caption>**Table 1**

Overview of All Conducted Runs: The Used Convective Mixing Schemes (EVD/EDMF/EDMF<sub>α</sub>/EDMF<sub>ρ</sub>/EDMF<sub>e</sub> (Details Can Be Found in Section 2.1.1)) and the Used Surface Freshwater Forcing</caption><thead><tr><th></th><th>2007–2011</th><th>2012–2017</th></tr></thead><tbody><tr><td>Convective mixing schemes (Section 3.1)</td><td>EVD, EDMF, EDMF<sub>α</sub>, EDMF<sub>ρ</sub>, EDMF<sub>e</sub> 100% precipitation</td><td></td></tr><tr><td rowspan="5">Perturbed freshwater forcing experiments (Section 3.2.1)</td><td rowspan="5">EDMF<sub>e</sub> 100% precipitation</td><td>EDMF<sub>e</sub>–40% precipitation</td></tr><tr><td>EDMF<sub>e</sub>–30% precipitation</td></tr><tr><td>EDMF<sub>e</sub>–20% precipitation</td></tr><tr><td>EDMF<sub>e</sub>–10% precipitation</td></tr><tr><td>EDMF<sub>e</sub>–0% precipitation</td></tr><tr><td rowspan="7">Convective mixing schemes combined with freshwater forcing experiments (Section 3.2.3)</td><td rowspan="7">EDMF<sub>e</sub> 100% precipitation</td><td>EVD–20% precipitation</td></tr><tr><td>EDMF–20% precipitation</td></tr><tr><td>EDMF<sub>α</sub>–20% precipitation</td></tr><tr><td>EDMF<sub>ρ</sub>–20% precipitation</td></tr><tr><td>EDMF<sub>e</sub>–20% precipitation</td></tr><tr><td>EVD–10% precipitation</td></tr><tr><td>EVD–0% precipitation</td></tr></tbody></table>

are carried out to investigate if a non-penetrative convection scheme with less freshwater can lead to similar convection patterns as a penetrative convection scheme.

## 2.2. Validation Using Observational Data

We used EN4.2.2 profiles (Good et al., 2013) spanning the period 1980–2022 and CTD data sets to create a monthly climatology for the region 66°S–64°S and 2°W–8°E (as used in Wilson et al. (2019)) in the upper 600 m (see Figure 1, green box). Additional CTD data from the following cruises were also included: MaudNESS (2005, (De Steur et al., 2007)) and SO-CHIC (2022, (Gülk et al., 2023)). For the climatology, EN4.2.2 profiles shallower than 300 m or with a salinity <32 g kg<sup>-1</sup> or temperature >10°C are excluded. The remaining profiles are vertically binned to 5 m intervals, and the MLD is defined by the density threshold Δσ = 0.01 kg m<sup>-3</sup> with respect to 10 m depth, as has been used by for example, Heuzé et al. (2015).

To compare the SIC of the simulations to observations, we use satellite data from AMSR-E and AMSR2 (Melsheimer & Sreen, 2019; Sreen et al., 2008) with 6.25 km horizontal resolution.

## 2.3. Model Analysis

For the vertical convective mixing scheme runs (Section 2.1.1), we estimate the ice import and export from the forced boundaries into the model domain. The ice import/export is crucial for the domain's freshwater balance. Therefore, we compute the ice transport $T_{ice}$ across the eastern, western and northern boundaries at the transition from the forced boundaries to the inner model domain. The southern boundary is confined by land, and there is no external ice forcing. The ice transport $T_{ice}$ is computed as

$$T_{ice} = \int_{bdy} h_{ice} C_{ice} \vec{u}_{ice} \cdot d\vec{n}, \qquad (5)$$

with the ice thickness $h_{ice}$, ice concentration $C_{ice}$ and ice velocity $\vec{u}_{ice}$ integrated across each boundary. Here, $d\vec{n}$ is oriented perpendicular to the boundary. The ice transport is converted into an equivalent freshening $\Delta S_{ice}$ for the upper ocean:

$$\Delta S_{ice} = \frac{\Delta S_{i} \rho_{ice} \int \Delta T_{ice} dt}{AH \rho_{w}}, \qquad (6)$$
<table><caption>**Table 1**

Overview of All Conducted Runs: The Used Convective Mixing Schemes (EVD/EDMF/EDMF<sub>α</sub>/EDMF<sub>ρ</sub>/EDMF<sub>e</sub> (Details Can Be Found in Section 2.1.1)) and the Used Surface Freshwater Forcing</caption><thead><tr><th></th><th>2007–2011</th><th>2012–2017</th></tr></thead><tbody><tr><td>Convective mixing schemes (Section 3.1)</td><td>EVD, EDMF, EDMF<sub>α</sub>, EDMF<sub>ρ</sub>, EDMF<sub>e</sub> 100% precipitation</td><td></td></tr><tr><td rowspan="5">Perturbed freshwater forcing experiments (Section 3.2.1)</td><td rowspan="5">EDMF<sub>e</sub> 100% precipitation</td><td>EDMF<sub>e</sub>–40% precipitation</td></tr><tr><td>EDMF<sub>e</sub>–30% precipitation</td></tr><tr><td>EDMF<sub>e</sub>–20% precipitation</td></tr><tr><td>EDMF<sub>e</sub>–10% precipitation</td></tr><tr><td>EDMF<sub>e</sub>–0% precipitation</td></tr><tr><td rowspan="7">Convective mixing schemes combined with freshwater forcing experiments (Section 3.2.3)</td><td rowspan="7">EDMF<sub>e</sub> 100% precipitation</td><td>EVD–20% precipitation</td></tr><tr><td>EDMF–20% precipitation</td></tr><tr><td>EDMF<sub>α</sub>–20% precipitation</td></tr><tr><td>EDMF<sub>ρ</sub>–20% precipitation</td></tr><tr><td>EDMF<sub>e</sub>–20% precipitation</td></tr><tr><td>EVD–10% precipitation</td></tr><tr><td>EVD–0% precipitation</td></tr></tbody></table>

are carried out to investigate if a non-penetrative convection scheme with less freshwater can lead to similar convection patterns as a penetrative convection scheme.

## 2.2. Validation Using Observational Data

We used EN4.2.2 profiles (Good et al., 2013) spanning the period 1980–2022 and CTD data sets to create a monthly climatology for the region 66°S–64°S and 2°W–8°E (as used in Wilson et al. (2019)) in the upper 600 m (see Figure 1, green box). Additional CTD data from the following cruises were also included: MaudNESS (2005, De Steur et al., 2007)) and SO-CHIC (2022, (Gülk et al., 2023)). For the climatology, EN4.2.2 profiles shallower than 300 m or with a salinity <32 g kg<sup>-1</sup> or temperature >10°C are excluded. The remaining profiles are vertically binned to 5 m intervals, and the MLD is defined by the density threshold Δσ = 0.01 kg m<sup>-3</sup> with respect to 10 m depth, as has been used by for example, Heuzé et al. (2015).

To compare the SIC of the simulations to observations, we use satellite data from AMSR-E and AMSR2 (Melsheimer & Sreen, 2019; Sreen et al., 2008) with 6.25 km horizontal resolution.

## 2.3. Model Analysis

For the vertical convective mixing scheme runs (Section 2.1.1), we estimate the ice import and export from the forced boundaries into the model domain. The ice import/export is crucial for the domain's freshwater balance. Therefore, we compute the ice transport $T_{ice}$ across the eastern, western and northern boundaries at the transition from the forced boundaries to the inner model domain. The southern boundary is confined by land, and there is no external ice forcing. The ice transport $T_{ice}$ is computed as

$$T_{ice} = \int_{bdy} h_{ice} C_{ice} \vec{u}_{ice} \cdot d\vec{n}, \qquad (5)$$

with the ice thickness $h_{ice}$, ice concentration $C_{ice}$ and ice velocity $\vec{u}_{ice}$ integrated across each boundary. Here, $d\vec{n}$ is oriented perpendicular to the boundary. The ice transport is converted into an equivalent freshening $\Delta S_{ice}$ for the upper ocean:

$$\Delta S_{ice} = \frac{\Delta S_{i} \rho_{ice} \int \Delta T_{ice} dt}{AH \rho_{w}}, \qquad (6)$$
where the difference in salinity between the ocean and sea ice is denoted by $\Delta S_i = 30 \text{ kg kg}^{-1}$, the density of sea ice as $\rho_{ice} = 916 \text{ kg m}^{-3}$, the area of the inner model domain as A, the depth of the upper ocean is set to H = 50 m, and the density of seawater to $\rho_w = 1,026 \text{ kg m}^{-3}$. The choice of H refers to the mean MLD estimated from the model output, which is around 60 m.

For all model runs, we compute the trend in the freshwater content, which we calculate as

$$F_c = \int \left( \frac{S_{ref} - S}{S_{ref}} \right) dV, \qquad (7)$$

with $S_{ref}$ being a reference salinity of 35 $\text{kg kg}^{-1}$ and S the salinity of the water integrated over the volume V of the inner model domain. The resulting time series for the freshwater content is used to derive the freshwater content trend via a linear regression.

We define deep convection events with a MLD threshold of $\ge 830 \text{ m}$ depth. This threshold is chosen in accordance with the observation of traces of deep-reaching convective mixing at greater depths following the 2016 MRP (Campbell et al., 2019), with full mixing over at least the top 800 m of the water column in 2017 (Cheon & Gordon, 2019) and to align with the model's vertical levels.

To quantify the water column's stratification and its potential to convect within the upper 830 m, we estimate the salt deficit $SD*$ and the convective resistance $CR$. The salt deficit measures the sea-ice growth required to trigger deep convection (Martinson & Iannuzzi, 1998):

$$SD* = \frac{h_{mld}}{\Delta S_i} \left[ \Delta S - \frac{\alpha}{\beta} \Delta \theta \right], \qquad (8)$$

with $h_{mld}$ as the MLD, and $\Delta S_i$ as the difference between the salinities of the ocean and sea ice. $\Delta S = S_d - S_{ml}$, using the salinities at a reference depth, $S_d$, and in the mixed layer, $S_{ml}$. Here, $\Delta \theta = \theta_d - \theta_f$ is the difference between the temperature at the reference depth, $\theta_d$, and that at the freezing point, $\theta_f$. The values $\alpha = 0.43 \cdot 10^{-4} \text{ K}^{-1}$ and $\beta = 7 \cdot 10^{-4} \text{ kg g}^{-1}$ are taken constant. We use the deep convection threshold of 830 m as reference depth. (The Supporting Information provides the full derivation of Equation 8.) The convective resistance $CR$ measures the buoyancy loss required to overturn the water column to a given depth H. We use the same approach as followed by Campbell et al. (2019):

$$CR = \frac{g}{\rho_w} \int_{0}^{H} (\sigma_{0}(H) - \sigma_{0}(z)) dz. \qquad (9)$$

Here, we choose H = 830 m, the same threshold as used for deep convection and the salt deficit; $\sigma_0$ as potential density anomaly with respect to a reference pressure of 0 dbar; and g = 9.81 $\text{ m s}^{-2}$ as gravitational acceleration due to gravity.

We investigate local differences between the Taylor Cap and the Maud Rise Flanks. We define the Taylor Cap to be located over the center of Maud Rise, where the bathymetry is shallower than 2,500 m. The Flanks are defined to encircle Maud Rise fully, therefore the 3,500 m isobath is smoothed and expanded by ~55 km. In the north of Maud Rise, the Flanks are limited by the 5,000 m isobath. This definition encircles the area of the warm water Halo (Figure 1, blue and orange encircled areas).

# 3. Results

## 3.1. The Vertical Convective Mixing Scheme

In this section, we analyze the results from the simulations employing EVD, EDMF and the modified versions of EDMF ($EDMF_\alpha$, EDMF$_\rho$, EDMF$_e$) for the period 2007–2017. The monthly climatological profiles of temperature, salinity and MLD for the region 66°S–64°S and 2°W–8°E (Figure 1, green box) are derived from these runs. Figure 2 depicts the mean annual cycles of temperature and salinity down to 600 m depth, comparing the simulations with EVD, EDMF and EDMF$_e$ to the observations. None of the runs display any trace of a polynya opening (see Figure S2 in Supporting Information S1). However, the runs do capture the general state of the
where the difference in salinity between the ocean and sea ice is denoted by $\Delta S_i = 30 \text{ g kg}^{-1}$, the density of sea ice as $\rho_{ice} = 916 \text{ kg m}^{-3}$, the area of the inner model domain as A, the depth of the upper ocean is set to H = 50 m, and the density of seawater to $\rho_w = 1,026 \text{ kg m}^{-3}$. The choice of H refers to the mean MLD estimated from the model output, which is around 60 m.

For all model runs, we compute the trend in the freshwater content, which we calculate as

$$F_c = \int \left( \frac{S_{ref} - S}{S_{ref}} \right) dV, \qquad (7)$$

with $S_{ref}$ being a reference salinity of 35 g kg$^{-1}$ and S the salinity of the water integrated over the volume V of the inner model domain. The resulting time series for the freshwater content is used to derive the freshwater content trend via a linear regression.

We define deep convection events with a MLD threshold of ≥830 m depth. This threshold is chosen in accordance with the observation of traces of deep-reaching convective mixing at greater depths following the 2016 MRP (Campbell et al., 2019), with full mixing over at least the top 800 m of the water column in 2017 (Cheon & Gordon, 2019) and to align with the model's vertical levels.

To quantify the water column's stratification and its potential to convect within the upper 830 m, we estimate the salt deficit SD* and the convective resistance CR. The salt deficit measures the sea-ice growth required to trigger deep convection (Martinson & Iannuzzi, 1998):

$$SD* = \frac{h_{mld}}{\Delta S_i} \left[ \Delta S - \frac{\alpha}{\beta} \Delta \theta \right], \qquad (8)$$

with $h_{mld}$ as the MLD, and $\Delta S_i$ as the difference between the salinities of the ocean and sea ice. $\Delta S = S_d - S_{ml}$, using the salinities at a reference depth, $S_d$, and in the mixed layer, $S_{ml}$. Here, $\Delta \theta = \theta_d - \theta_f$ is the difference between the temperature at the reference depth, $\theta_d$, and that at the freezing point, $\theta_f$. The values $\alpha = 0.43 \cdot 10^{-4} \text{ K}^{-1}$ and $\beta = 7 \cdot 10^{-4} \text{ kg g}^{-1}$ are taken constant. We use the deep convection threshold of 830 m as reference depth. (The Supporting Information provides the full derivation of Equation 8.) The convective resistance CR measures the buoyancy loss required to overturn the water column to a given depth H. We use the same approach as followed by Campbell et al. (2019):

$$CR = \frac{g}{\rho_w} \int_{0}^{H} (\sigma_{0}(H) - \sigma_{0}(z)) dz. \qquad (9)$$

Here, we choose H = 830 m, the same threshold as used for deep convection and the salt deficit; $\sigma_0$ as potential density anomaly with respect to a reference pressure of 0 dbar; and g = 9.81 m s$^{-2}$ as gravitational acceleration due to gravity.

We investigate local differences between the Taylor Cap and the Maud Rise Flanks. We define the Taylor Cap to be located over the center of Maud Rise, where the bathymetry is shallower than 2,500 m. The Flanks are defined to encircle Maud Rise fully, therefore the 3,500 m isobath is smoothed and expanded by ~55 km. In the north of Maud Rise, the Flanks are limited by the 5,000 m isobath. This definition encircles the area of the warm water Halo (Figure 1, blue and orange encircled areas).

# 3. Results

## 3.1. The Vertical Convective Mixing Scheme

In this section, we analyze the results from the simulations employing EVD, EDMF and the modified versions of EDMF ($EDMF_\alpha$, EDMF$_\rho$, EDMF$_e$) for the period 2007–2017. The monthly climatological profiles of temperature, salinity and MLD for the region 66°S–64°S and 2°W–8°E (Figure 1, green box) are derived from these runs. Figure 2 depicts the mean annual cycles of temperature and salinity down to 600 m depth, comparing the simulations with EVD, EDMF and EDMF$_e$ to the observations. None of the runs display any trace of a polynya opening (see Figure S2 in Supporting Information S1). However, the runs do capture the general state of the
<figure><img src="image_5.png" /><figcaption>**Figure 2.** (a) Temperature and (b) salinity of the upper 600 m in the region 66°S–64°S; 2°W–8°E from observations, and difference between the simulations and observations (c–h). Red dashed (black solid) lines indicate the mixed layer depth from the simulations (observations).</figcaption></figure>

region for non-polynya years. The observations are also largely obtained during non-polynya years, such that a comparison is reasonable. Comparing temperature fields shows that a cold bias occurs in the upper ocean, and a warm bias below the MLD, from January to April in the simulations. During July to November a warm bias occurs at the base of the MLD in the simulation using EVD. This warm bias vanishes when using EDMF, and is reduced when using $EDMF_c$. These characteristics can also be found in the root-mean-square-error (RMSE) of temperature for the upper 640 m, where EVD leads to a RMSE of 0.27°C, EDMF to 0.33°C, $EDMF_α$ to 0.33°C, $EDMF_ρ$ to 0.34°C and $EDMF_c$ to 0.30°C. The simulation-based climatologies also exhibit and overly fresh surface layer throughout the year, and a fresh bias below the MLD in February to July. Here, the RMSE of salinity is for all simulations around $0.16 \text{ g kg}^{-1}$. The representation of the MLD varies between the seasons within the schemes. From January to July, EVD shows smaller differences compared to observation (Figures 2c and 2d), while in August to November, the polynya season, $EDMF_c$ performs better. Thus, the smallest RMSE for the MLD is found when using $EDMF_c$ with 10.92 m, and the largest RMSE is found for $EDMF_ρ$ with 16.61 m. The RMSE of the MLD for EVD is 13.3 m, for EDMF is 15.04 m, and for $EDMF_α$ is 14.63 m. Considering the temperature and salinity biases, this indicates that the upper layer is too cold and fresh in all simulations, and thus that the modeled upper ocean is too stratified.

To quantify the origin of the temperature and salinity biases, the equivalent freshening $\Delta S_{ice}$ (Equation 6) and the trend of the freshwater content $F_c$ (Equation 7) are estimated for each simulation (Table 2). The lowest freshening
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) Temperature and (b) salinity of the upper 600 m in the region 66°S–64°S; 2°W–8°E from observations, and difference between the simulations and observations (c–h). Red dashed (black solid) lines indicate the mixed layer depth from the simulations (observations).</figcaption></figure>

region for non-polynya years. The observations are also largely obtained during non-polynya years, such that a comparison is reasonable. Comparing temperature fields shows that a cold bias occurs in the upper ocean, and a warm bias below the MLD, from January to April in the simulations. During July to November a warm bias occurs at the base of the MLD in the simulation using EVD. This warm bias vanishes when using EDMF, and is reduced when using $EDMF_c$. These characteristics can also be found in the root-mean-square-error (RMSE) of temperature for the upper 640 m, where EVD leads to a RMSE of 0.27°C, EDMF to 0.33°C, $EDMF_\alpha$ to 0.33°C, $EDMF_\rho$ to 0.34°C and $EDMF_e$ to 0.30°C. The simulation-based climatologies also exhibit and overly fresh surface layer throughout the year, and a fresh bias below the MLD in February to July. Here, the RMSE of salinity is for all simulations around $0.16 \text{ g kg}^{-1}$. The representation of the MLD varies between the seasons within the schemes. From January to July, EVD shows smaller differences compared to observation (Figures 2c and 2d), while in August to November, the polynya season, $EDMF_e$ performs better. Thus, the smallest RMSE for the MLD is found when using $EDMF_e$ with 10.92 m, and the largest RMSE is found for $EDMF_\rho$ with 16.61 m. The RMSE of the MLD for EVD is 13.3 m, for EDMF is 15.04 m, and for $EDMF_\alpha$ is 14.63 m. Considering the temperature and salinity biases, this indicates that the upper layer is too cold and fresh in all simulations, and thus that the modeled upper ocean is too stratified.

To quantify the origin of the temperature and salinity biases, the equivalent freshening $\Delta S_{ice}$ (Equation 6) and the trend of the freshwater content $F_c$ (Equation 7) are estimated for each simulation (Table 2). The lowest freshening
<table><caption>**Table 2**
Equivalent Freshening $\Delta S_{ice}$, and Freshwater Content Trend, for the Model Runs Over the Period 2007 to 2017</caption><thead><tr><th></th><th>ΔS<sub>ice</sub> [g kg<sup>-1</sup>]</th><th>Trend of F<sub>c</sub> [10<sup>-3</sup> Sv]</th></tr></thead><tbody><tr><td>EVD</td><td>-1.00</td><td>1.19 ± 0.06</td></tr><tr><td>EDMF</td><td>-2.58</td><td>1.37 ± 0.06</td></tr><tr><td>EDMF<sub>α</sub></td><td>-2.55</td><td>1.36 ± 0.06</td></tr><tr><td>EDMF<sub>ρ</sub></td><td>-2.86</td><td>1.18 ± 0.06</td></tr><tr><td>EDMF<sub>c</sub></td><td>-2.16</td><td>1.20 ± 0.06</td></tr></tbody></table>

is found when using EVD, and is highest in the simulation with EDMF<sub>c</sub>. In contrast, EDMF<sub>ρ</sub> shows the lowest freshwater content trend, which is slightly lower than in the runs using EVD and EDMF<sub>c</sub>, even though the freshening rates of the EDMF simulations are considerably higher than in the EVD run. These differences indicate that the model domain accumulates ice, which then leads to an overly stratified upper ocean. The accumulation of ice in the domain results partly from the open-boundary conditions of the configuration. The choice of the convective mixing scheme does not lead to a major change in the ocean circulation or sea-ice motion. A winter-mean sea-ice vorticity and a meridional cross section of zonal velocity can be found in Figure S3 in Supporting Information S1. A region with minor differences with increased zonal velocities and a more negative sea-ice vorticity occurs north of the

Halo, but these changes cannot explain the observed differences in $\Delta S_{ice}$. Therefore, the increased $\Delta S_{ice}$ is likely related to the reduced heat bias found in the climatologies (Figures 2e and 2g). Comparing across all runs: the MLD during the polynya season, the smallest RMSE for the MLD, comparable RMSE values for temperature and salinity, and the trend in the freshwater content are best represented using EDMF<sub>c</sub>. Thus, we employ EDMF<sub>c</sub> for performing the sensitivity runs.

## 3.2. Sensitivity Runs

From the five runs performed with reduced precipitation and using EDMF<sub>c</sub> (Section 2.1.2), we identify the extent of reduction necessary to best reproduce a polynya opening at Maud Rise. Accordingly, we compute the trend in freshwater content and the number of days with a SIC <20% between 1st July and 30th September of each year (Section 3.2.1). Based on these results, the run with strongest agreement with observations is selected and used to investigate the processes leading to deep convection and polynya opening (Section 3.2.2).

### 3.2.1. The Impact of Surface Freshwater Forcing

The trend of the freshwater content in all runs using perturbed surface freshwater forcing including the spin-up is $\approx 0.9 \cdot 10^{-3}$ Sv (the exact values can be found in the Table S1 in Supporting Information S1). This trend is $\approx 0.3 \cdot 10^{-3}$ Sv lower than in the EDMF<sub>c</sub> simulation in the first set of tests (Table 2). The reduction of precipitation also leads to a modification in the ice import and export, as well as the ice production, such that the previously found ice accumulation in the domain is now compensated for.

The number of days with a SIC <20% between 1st July and 30th September of each year are estimated for each run (Figure 3). According to previous observational studies of the region, a brief polynya is expected in 2016 and a longer one in 2017 at Maud Rise, with no major openings in the period 2012–2015. A minor opening of 2–3 days' duration occurred in 2013 (Campbell et al., 2019). The polynya openings of 2016 and 2017 are found at the northeastern flank of Maud Rise in all runs using 30% precipitation or less. The simulation using 40% precipitation does not show a signature of a sea-ice opening in the region. The simulations with 10% precipitation or less display polynya openings even earlier than 2016. The 20% simulation exhibits a polynya in 2016 and 2017 at Maud Rise, with a weak opening in 2013. The signature is stronger in 2017 than in 2016, which agrees with expectations from observations. The polynya locations are not exactly those found in the observations, but are in the polynya-prone region. On the basis of the trend of the freshwater content and realistic SIC patterns in the 20% precipitation simulation, we select that simulation to investigate the mechanisms controlling polynya formation and deep convection next.

### 3.2.2. Polynya Formation and Deep Convection Events

The 20% precipitation simulation is analyzed to elucidate the processes governing polynya formation and deep convection events. First, the polynya openings in 2016 and 2017 are compared to satellite observations. Then the locations and cause of deep convection are determined. In the observations, the polynya opens on the 25th July 2016 and on 3rd September 2017, after a full ice cover in the area. In the simulation, the openings occur in late June/early July in both years, and in regions that have only been covered by a very thin ice layer or are still open ocean (Figure 4). At this time the northeastern part of the domain is not ice-covered. The shapes of the 10% SIC contour (red contour) are broadly similar between observations and simulation, for example, on 6th August 2016
<table><caption>**Table 2**
Equivalent Freshening $\Delta S_{ice}$, and Freshwater Content Trend, for the Model Runs Over the Period 2007 to 2017</caption><thead><tr><th></th><th>ΔS<sub>ice</sub> [g kg<sup>-1</sup>]</th><th>Trend of F<sub>c</sub> [10<sup>-3</sup> Sv]</th></tr></thead><tbody><tr><td>EVD</td><td>-1.00</td><td>1.19 ± 0.06</td></tr><tr><td>EDMF</td><td>-2.58</td><td>1.37 ± 0.06</td></tr><tr><td>EDMF<sub>α</sub></td><td>-2.55</td><td>1.36 ± 0.06</td></tr><tr><td>EDMF<sub>ρ</sub></td><td>-2.86</td><td>1.18 ± 0.06</td></tr><tr><td>EDMF<sub>c</sub></td><td>-2.16</td><td>1.20 ± 0.06</td></tr></tbody></table>

is found when using EVD, and is highest in the simulation with EDMF<sub>ρ</sub>. In contrast, EDMF<sub>ρ</sub> shows the lowest freshwater content trend, which is slightly lower than in the runs using EVD and EDMF<sub>c</sub>, even though the freshening rates of the EDMF simulations are considerably higher than in the EVD run. These differences indicate that the model domain accumulates ice, which then leads to an overly stratified upper ocean. The accumulation of ice in the domain results partly from the open-boundary conditions of the configuration. The choice of the convective mixing scheme does not lead to a major change in the ocean circulation or sea-ice motion. A winter-mean sea-ice vorticity and a meridional cross section of zonal velocity can be found in Figure S3 in Supporting Information S1. A region with minor differences with increased zonal velocities and a more negative sea-ice vorticity occurs north of the

Halo, but these changes cannot explain the observed differences in $\Delta S_{ice}$. Therefore, the increased $\Delta S_{ice}$ is likely related to the reduced heat bias found in the climatologies (Figures 2e and 2g). Comparing across all runs: the MLD during the polynya season, the smallest RMSE for the MLD, comparable RMSE values for temperature and salinity, and the trend in the freshwater content are best represented using EDMF<sub>c</sub>. Thus, we employ EDMF<sub>c</sub> for performing the sensitivity runs.

## 3.2. Sensitivity Runs

From the five runs performed with reduced precipitation and using EDMF<sub>c</sub> (Section 2.1.2), we identify the extent of reduction necessary to best reproduce a polynya opening at Maud Rise. Accordingly, we compute the trend in freshwater content and the number of days with a SIC <20% between 1st July and 30th September of each year (Section 3.2.1). Based on these results, the run with strongest agreement with observations is selected and used to investigate the processes leading to deep convection and polynya opening (Section 3.2.2).

### 3.2.1. The Impact of Surface Freshwater Forcing

The trend of the freshwater content in all runs using perturbed surface freshwater forcing including the spin-up is $\approx 0.9 \cdot 10^{-3}$ Sv (the exact values can be found in the Table S1 in Supporting Information S1). This trend is $\approx 0.3 \cdot 10^{-3}$ Sv lower than in the EDMF<sub>c</sub> simulation in the first set of tests (Table 2). The reduction of precipitation also leads to a modification in the ice import and export, as well as the ice production, such that the previously found ice accumulation in the domain is now compensated for.

The number of days with a SIC <20% between 1st July and 30th September of each year are estimated for each run (Figure 3). According to previous observational studies of the region, a brief polynya is expected in 2016 and a longer one in 2017 at Maud Rise, with no major openings in the period 2012–2015. A minor opening of 2–3 days' duration occurred in 2013 (Campbell et al., 2019). The polynya openings of 2016 and 2017 are found at the northeastern flank of Maud Rise in all runs using 30% precipitation or less. The simulation using 40% precipitation does not show a signature of a sea-ice opening in the region. The simulations with 10% precipitation or less display polynya openings even earlier than 2016. The 20% simulation exhibits a polynya in 2016 and 2017 at Maud Rise, with a weak opening in 2013. The signature is stronger in 2017 than in 2016, which agrees with expectations from observations. The polynya locations are not exactly those found in the observations, but are in the polynya-prone region. On the basis of the trend of the freshwater content and realistic SIC patterns in the 20% precipitation simulation, we select that simulation to investigate the mechanisms controlling polynya formation and deep convection next.

### 3.2.2. Polynya Formation and Deep Convection Events

The 20% precipitation simulation is analyzed to elucidate the processes governing polynya formation and deep convection events. First, the polynya openings in 2016 and 2017 are compared to satellite observations. Then the locations and cause of deep convection are determined. In the observations, the polynya opens on the 25th July 2016 and on 3rd September 2017, after a full ice cover in the area. In the simulation, the openings occur in late June/early July in both years, and in regions that have only been covered by a very thin ice layer or are still open ocean (Figure 4). At this time the northeastern part of the domain is not ice-covered. The shapes of the 10% SIC contour (red contour) are broadly similar between observations and simulation, for example, on 6th August 2016
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Number of days with SIC <20% between 1st July and 30th September of each year, per sensitivity run and per year. The gray lines indicate the 2,500 m, 3,500 m, 5,000 m isobaths and the black contour encircles the area of the Taylor Cap.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Snapshots of SIC from satellite observations (a–d) and the model run (e–h) in 2016 and 2017. The 50% (10%) SIC contour is indicated in blue (red). (i) Area with SIC <50% (blue; <10%; red) in the region 4°W–10°E and 62°S–67°S from the model (solid) and the observations (dotted).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Number of days with SIC <20% between 1st July and 30th September of each year, per sensitivity run and per year. The gray lines indicate the 2,500 m, 3,500 m, 5,000 m isobaths and the black contour encircles the area of the Taylor Cap.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Snapshots of SIC from satellite observations (a–d) and the model run (e–h) in 2016 and 2017. The 50% (10%) SIC contour is indicated in blue (red). (i) Area with SIC <50% (blue; <10%; red) in the region 4°W–10°E and 62°S–67°S from the model (solid) and the observations (dotted).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Locations of initial deep convection (yellow triangles) in panel (a) 2016 and (b) 2017 from the simulation. The encircled area shows the 50% SIC contour from satellite data (red lines) and simulation (blue lines) at the first occurrence of the polynya. The background shows the bathymetry with the Taylor Cap encircled (black contour).</figcaption></figure>

(cf. Figures 4b and 4f) and on 25th September 2017 (observations, Figure 4d) and 12th July 2017 (model, Figure 4h). The location of this contour indicates Maud Rise as an anchor point for the polynya. The 50% SIC contour (blue line) shows that in the simulation the polynya is present over a larger area. This can be seen in the time series of area below SIC <50% in the region 4°W–10°E and 62°S–67°S (Figure 4i). This area is chosen to exclude the delayed ice formation east of 10°E in the model, which appears in some of the winter seasons (e.g., Figure 3), and open-water areas due coastal polynyas south of 67°S. This time series shows that the simulation and observations of the area with SIC <10% converge in August 2016. This is also true for the area with SIC <50%, but following the polynya event the polynya stays open in the simulation while it closes down in the observations. In 2017, the modeled polynya appears too early, with a large area until September 2017, while in the observations the opening appears from September onward.

The locations of deep convection in 2016 and 2017 are identified in the simulation by using the MLD threshold ≥830 m. The initial deep convection locations are compared to the position of the polynya opening in the observations and the model by using the 50% SIC contour (Figure 5). The comparison of deep convection sites in the model to satellite sea-ice openings is reasonable, since polynyas and deep convection are often concurrent and physically linked. The limited spatio-temporal density of observations does not allow for an estimate of the onset of deep convection or region in 2016 and 2017. In 2016, the simulation shows the first signals of deep convection on the 17th July at the northern flank of Maud Rise, between the Taylor Cap and the Halo. The satellite data reveals an opening of the sea ice on the 25th July, also north of Maud Rise (Figure 5a). The opening in the satellite data is 8 days later than the onset of deep convection in the simulation. At the onset of deep convection in the simulation, the Halo as well as the northeastern part of the model domain display a SIC below 50%, and the ice free area is larger than in the observations. In 2017, the deep convection occurs on the 30th June and 1st July. On the 30th June the deep convection signal is found in a large region on the northeastern flank of Maud Rise, and a day later at the western flank. Satellite data displays the opening on the 3rd September, slightly west of the locations of deep convection in the simulation (Figure 5b). Due to the early onset of deep convection, large parts of the northern model domain have a SIC below 50%. The offset in the timing of the opening may be related to some extent to the uncertainties in atmospheric reanalysis products, due to the limitations in observations as well as resolution of the used product.

The characteristics of deep convection across Maud Rise are investigated via a three-step approach. First, we determine the day of onset of deep convection at each grid cell in each year. Second, we estimate the time lag (relative to the onset of convection in the region) as the temporal difference between the local onset of convection and the first onset within the model domain (Figures 6a and 6b). And third, we estimate the duration of convection as the number of days in a year during which the MLD was deeper than the deep convection threshold of 830 m (Figures 6c and 6d). In 2016, initially only a few grid cells display deep convection. Subsequently deep convection
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Locations of initial deep convection (yellow triangles) in panel (a) 2016 and (b) 2017 from the simulation. The encircled area shows the 50% SIC contour from satellite data (red lines) and simulation (blue lines) at the first occurrence of the polynya. The background shows the bathymetry with the Taylor Cap encircled (black contour).</figcaption></figure>

(cf. Figures 4b and 4f) and on 25th September 2017 (observations, Figure 4d) and 12th July 2017 (model, Figure 4h). The location of this contour indicates Maud Rise as an anchor point for the polynya. The 50% SIC contour (blue line) shows that in the simulation the polynya is present over a larger area. This can be seen in the time series of area below SIC <50% in the region 4°W–10°E and 62°S–67°S (Figure 4i). This area is chosen to exclude the delayed ice formation east of 10°E in the model, which appears in some of the winter seasons (e.g., Figure 3), and open-water areas due coastal polynyas south of 67°S. This time series shows that the simulation and observations of the area with SIC <10% converge in August 2016. This is also true for the area with SIC <50%, but following the polynya event the polynya stays open in the simulation while it closes down in the observations. In 2017, the modeled polynya appears too early, with a large area until September 2017, while in the observations the opening appears from September onward.

The locations of deep convection in 2016 and 2017 are identified in the simulation by using the MLD threshold ≥830 m. The initial deep convection locations are compared to the position of the polynya opening in the observations and the model by using the 50% SIC contour (Figure 5). The comparison of deep convection sites in the model to satellite sea-ice openings is reasonable, since polynyas and deep convection are often concurrent and physically linked. The limited spatio-temporal density of observations does not allow for an estimate of the onset of deep convection or region in 2016 and 2017. In 2016, the simulation shows the first signals of deep convection on the 17th July at the northern flank of Maud Rise, between the Taylor Cap and the Halo. The satellite data reveals an opening of the sea ice on the 25th July, also north of Maud Rise (Figure 5a). The opening in the satellite data is 8 days later than the onset of deep convection in the simulation. At the onset of deep convection in the simulation, the Halo as well as the northeastern part of the model domain display a SIC below 50%, and the ice free area is larger than in the observations. In 2017, the deep convection occurs on the 30th June and 1st July. On the 30th June the deep convection signal is found in a large region on the northeastern flank of Maud Rise, and a day later at the western flank. Satellite data displays the opening on the 3rd September, slightly west of the locations of deep convection in the simulation (Figure 5b). Due to the early onset of deep convection, large parts of the northern model domain have a SIC below 50%. The offset in the timing of the opening may be related to some extent to the uncertainties in atmospheric reanalysis products, due to the limitations in observations as well as resolution of the used product.

The characteristics of deep convection across Maud Rise are investigated via a three-step approach. First, we determine the day of onset of deep convection at each grid cell in each year. Second, we estimate the time lag (relative to the onset of convection in the region) as the temporal difference between the local onset of convection and the first onset within the model domain (Figures 6a and 6b). And third, we estimate the duration of convection as the number of days in a year during which the MLD was deeper than the deep convection threshold of 830 m (Figures 6c and 6d). In 2016, initially only a few grid cells display deep convection. Subsequently deep convection
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Time lag between the initial convection locations (cyan areas) and the onset of deep convection in the respective grid cell (a, b) and duration (c, d) of deep convection in 2016 (left) and 2017 (right). The hatched area in panel (b) shows the extent of deep convection in 2016. Gray lines indicate the 2,500, 3,500 and 5,000 m isobaths and the black contour encircles the area of the Taylor Cap.</figcaption></figure>

slowly spreads across Maud Rise (Figure 6a). The time lag shows that, from the initial deep convection location, the signal then spreads along the direction of the retroflection current, as far as the ridge east of Maud Rise. This ridge is reached within a time lag of 10 days, which corresponds well with the advective time scale of the retroflection current. From the ridge, deep convection signals extend across the Taylor Cap and to the west of Maud Rise. The area of initial convection convected for a brief period (<5 days) only, whereas the region to the north of the ridge exhibits convection for up to 47 days. This agrees with the location of the polynya opening (Figure 4). In 2017, a large area of the Taylor Cap displays convection within 10 days of the first deep convection event. After such initial rapid expansion, deep convection continues to spread slowly, and is mainly confined to areas of the Taylor Cap and Halo. This spatial limitation in deep convection expansion may be related to the high horizontal resolution and adequate representation of the bathymetry, which might be too coarse in lower-resolution models and lead to excessively large polynya patterns. The rapid expansion of deep convection in 2017 is due to preconditioning by the deep convection event in 2016, which reduces ambient stratification across the region.

The above overview of deep convection in 2016 and 2017 reveals differences in the spatial extents of convection onset, as well as in the duration and time lag of convective events. This indicates that different processes could have triggered convection in each year. Thus, we next focus on the 30-day period preceding the onset of deep convection and examine the possible role of a suite of processes (e.g., thermobaricity, sea-ice growth, buoyancy loss and surface salt anomalies) in triggering convection. Instabilities within the water column are identified by the maximum convective plume velocity, $w_p$, in the region of the onset of deep convection (Figures 7a and 7b). The onset region is that indicated by a purple box in Figures 7c and 7d. In 2016, the maximum plume velocity is mainly found within the mixed layer and at depths of between 200 and 250 m prior to the initial deep convection event. This depth range approximately aligns with the interface between the WDW and the cold surface water. In 2017, the maximum plume velocity is mainly found within the mixed layer, below which only weak signals are found in the period preceding deep convection. The stratification of the water column and its potential to overturn can be assessed via the salt deficit $\text{SD}^*$(Equation 8, Figures 7c and 7d). The $\text{SD}^*$ is computed from daily output and averaged over the 30 days preceding the deep convection. This reveals that, in 2016, the stratification was
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Time lag between the initial convection locations (cyan areas) and the onset of deep convection in the respective grid cell (a, b) and duration (c, d) of deep convection in 2016 (left) and 2017 (right). The hatched area in panel (b) shows the extent of deep convection in 2016. Gray lines indicate the 2,500, 3,500 and 5,000 m isobaths and the black contour encircles the area of the Taylor Cap.</figcaption></figure>

slowly spreads across Maud Rise (Figure 6a). The time lag shows that, from the initial deep convection location, the signal then spreads along the direction of the retroflection current, as far as the ridge east of Maud Rise. This ridge is reached within a time lag of 10 days, which corresponds well with the advective time scale of the retroflection current. From the ridge, deep convection signals extend across the Taylor Cap and to the west of Maud Rise. The area of initial convection convected for a brief period (<5 days) only, whereas the region to the north of the ridge exhibits convection for up to 47 days. This agrees with the location of the polynya opening (Figure 4). In 2017, a large area of the Taylor Cap displays convection within 10 days of the first deep convection event. After such initial rapid expansion, deep convection continues to spread slowly, and is mainly confined to areas of the Taylor Cap and Halo. This spatial limitation in deep convection expansion may be related to the high horizontal resolution and adequate representation of the bathymetry, which might be too coarse in lower-resolution models and lead to excessively large polynya patterns. The rapid expansion of deep convection in 2017 is due to preconditioning by the deep convection event in 2016, which reduces ambient stratification across the region.

The above overview of deep convection in 2016 and 2017 reveals differences in the spatial extents of convection onset, as well as in the duration and time lag of convective events. This indicates that different processes could have triggered convection in each year. Thus, we next focus on the 30-day period preceding the onset of deep convection and examine the possible role of a suite of processes (e.g., thermobaricity, sea-ice growth, buoyancy loss and surface salt anomalies) in triggering convection. Instabilities within the water column are identified by the maximum convective plume velocity, $w_p$, in the region of the onset of deep convection (Figures 7a and 7b). The onset region is that indicated by a purple box in Figures 7c and 7d. In 2016, the maximum plume velocity is mainly found within the mixed layer and at depths of between 200 and 250 m prior to the initial deep convection event. This depth range approximately aligns with the interface between the WDW and the cold surface water. In 2017, the maximum plume velocity is mainly found within the mixed layer, below which only weak signals are found in the period preceding deep convection. The stratification of the water column and its potential to overturn can be assessed via the salt deficit $\text{SD}^*$(Equation 8, Figures 7c and 7d). The $\text{SD}^*$ is computed from daily output and averaged over the 30 days preceding the deep convection. This reveals that, in 2016, the stratification was
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Maximum plume velocity in the region of initial deep convection in panel (a) 2016 and (b) 2017. The white (black) dotted line indicates the onset of deep convection (30 days before). The gray line denotes the maximum mixed layer depth in the region of initial deep convection. (c, d) Mean SD\* and (e, f) mean convective resistance of the 30 days prior to the onset of deep convection in 2016 and 2017. (g, h) Salt anomaly in the upper 100 m, of the 30 days before convection onset in 2016 and 2017. Purple boxes in panels (c, d) indicate the regions of initial deep convection. Gray lines indicate the 2,500, 3,500 and 5,000 m isobaths and the black contour encircles the area of the Taylor Cap.</figcaption></figure>

reduced and only small sea-ice growth would be necessary to overturn the water column at the northern flank of Maud Rise. In 2017, the mean SD\* is reduced compared to 2016, and shows that the water column has been preconditioned the previous year through an overall reduction of stratification. The low mean SD\* values are aligned with the area in which rapid overturning occurs after the onset of deep convection. The convective resistance CR (Equation 9, Figures 7e and 7f) assesses the buoyancy loss required to overturn the upper water column as a temporal mean of the 30 days preceding the deep convection. Here we use a reference depth of 830 m. The pattern resembles that of the SD\*, with a low resistance at the northern flank of Maud Rise in 2016 and an overall weaker resistance in the region of Maud Rise in 2017. This is especially so at locations that have been
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Maximum plume velocity in the region of initial deep convection in panel (a) 2016 and (b) 2017. The white (black) dotted line indicates the onset of deep convection (30 days before). The gray line denotes the maximum mixed layer depth in the region of initial deep convection. (c, d) Mean SD\* and (e, f) mean convective resistance of the 30 days prior to the onset of deep convection in 2016 and 2017. (g, h) Salt anomaly in the upper 100 m, of the 30 days before convection onset in 2016 and 2017. Purple boxes in panels (c, d) indicate the regions of initial deep convection. Gray lines indicate the 2,500, 3,500 and 5,000 m isobaths and the black contour encircles the area of the Taylor Cap.</figcaption></figure>

reduced and only small sea-ice growth would be necessary to overturn the water column at the northern flank of Maud Rise. In 2017, the mean SD\* is reduced compared to 2016, and shows that the water column has been preconditioned the previous year through an overall reduction of stratification. The low mean SD\* values are aligned with the area in which rapid overturning occurs after the onset of deep convection. The convective resistance CR (Equation 9, Figures 7e and 7f) assesses the buoyancy loss required to overturn the upper water column as a temporal mean of the 30 days preceding the deep convection. Here we use a reference depth of 830 m. The pattern resembles that of the SD\*, with a low resistance at the northern flank of Maud Rise in 2016 and an overall weaker resistance in the region of Maud Rise in 2017. This is especially so at locations that have been
preconditioned. The mean salt anomaly in the upper 100 m for the 30 days preceding convection is computed from the difference between the salinity at that time and the temporal mean salinity for the period 2016–2017, and then averaged over the upper 100 m (Figures 7g and 7h). This depth range is used following the method outlined by Kurtakoti et al. (2018). In 2016, the mean salinity anomaly is modest in the region of Maud Rise, with the largest anomalies found north of the Halo and on the Rise’s southwestern flank. This is contrary to the patterns in 2017, when large anomalies occur on the north and northeastern flanks, especially in the locations in which deep convection is initiated.

In summary, the deep convection in 2016 resulted from instabilities at depth in combination with a reduced stratification. In turn, convection in 2017 was preconditioned by the 2016 event, and triggered by an anomalous salt advection within the upper 100 m.

### 3.2.3. Impact of the Vertical Convective Mixing Scheme

To understand the role of the vertical convective mixing scheme adopted in the model in shaping polynya formation and convective activity, we repeat the sensitivity run with 20% precipitation using the different convective mixing schemes (EVD, EDMF and $EDMF_ρ$) after the spin-up period. Since the vertical convective mixing scheme is the only change in the model setup, any differences must relate to the choice of scheme. We focus on the differences in the maximum MLD, the stratification below 50 m beneath the MLD (hereafter $MLD_{+50}$), and relation of $N^2$ to the velocity of the convective plumes in 2016. We use the $MLD_{+50}$ threshold, to exclude mixing within the often highly turbulent mixed layer.

In all runs, a maximum MLD deeper than 800 m is found in 2016 (Figures 8a–8d), but the expansiveness of the area with such deep mixed layers differs between simulations. Using EVD leads to deep MLDs at the northern flank of Maud Rise, as well as toward the south (at 4°W, 68.5°S) in small, confined regions (Figure 8a). The maximum MLD patterns are identical between EDMF and $EDMF_α$ (not shown). Changing the thermal expansion and haline contraction coefficients to non-fixed values within the scheme has no impact on the deep MLD properties. Using $EDMF_ρ$ (Figure 8c) produces the largest region of deep MLDs amongst all runs. This results from changing the potential density to in situ density within the scheme. Including the frictional term (Equation 3), which is related to the lateral entrainment (EDMF$_e$, Figure 8d) leads to a smaller region of deep MLDs. The lateral entrainment term (Equation 4) is rather dominated by the turbulence at the plume-environment interface than the background entrainment $ε_g$. The turbulence term is on average for 2016 around 0.86, while the background entrainment is kept constant with 0.001. Further, it is notable that all EDMF schemes produce a clearer maximum MLD patterns, while with EVD the region south of Maud Rise (which host MLDs shallower than 200 m in all EDMF runs) exhibits MLDs of about 500 m. This difference in convective behavior likely stems from the distinct ways in which the EVD and EDMF schemes represent convective mixing. EDMF works on a sub-grid scale and penetrates through the water column, while EVD mixes neighboring grid cells vertically. In the runs using EVD and 10% and 0% precipitation (Figures 8e and 8f), the maximum MLD reaches below 800 m depth, localized on top of Maud Rise and its flanks. Additionally, deep MLDs (>500 m) can be found north of Astrid Ridge and south of Maud Rise. The reduction from 10% to 0% precipitation leads to an increased area with MLD >800 m and an overall increase in the MLD, but with consistent spatial patterns. To conclude, all tested EDMF schemes are able to produce localized patterns of deep convection that only impact specific regions, while the commonly used EVD generates weaker convection affecting a larger region.

The stratification in the deeper water column is also impacted by the scheme choice. (Monthly means of 2016 for temperature and salinity for the different convective mixing schemes can be found in the Figure S4 in Supporting Information S1.) Therefore, the seasonal-horizontal mean stratification for the Taylor Cap and Flanks (Figure 1 blue and orange encircled area) below $MLD_{+50}$ is estimated for each depth level from the daily output for the four seasons in 2016: summer (JFM), autumn (AMJ), winter (JAS) and spring (OND) (Figure 9). In summer, the distribution of $N^2$ below $MLD_{+50}$ in the Taylor Cap exhibits the largest range of $N^2$ values when using EVD. The EDMF-produced distribution resembles that for EVD, but the maximum (minimum) $N^2$ is smaller (larger). The smallest values of $N^2$ are found when using $EDMF_e$. The distribution of stratification at the Flanks is almost identical between the schemes. The median for both regions is just above 0. In autumn, the largest $N^2$ values for the Taylor Cap are found when using EDMF, and smallest for $EDMF_ρ$. The same patterns are found at the Flanks. The weak stratification at the Flanks and Taylor Cap in $EDMF_ρ$ is likely affected by the onset of deep convection in late June, when the median is above 0. The stratification of the Taylor Cap is affected by the intrusions of colder
preconditioned. The mean salt anomaly in the upper 100 m for the 30 days preceding convection is computed from the difference between the salinity at that time and the temporal mean salinity for the period 2016–2017, and then averaged over the upper 100 m (Figures 7g and 7h). This depth range is used following the method outlined by Kurtakoti et al. (2018). In 2016, the mean salinity anomaly is modest in the region of Maud Rise, with the largest anomalies found north of the Halo and on the Rise’s southwestern flank. This is contrary to the patterns in 2017, when large anomalies occur on the north and northeastern flanks, especially in the locations in which deep convection is initiated.

In summary, the deep convection in 2016 resulted from instabilities at depth in combination with a reduced stratification. In turn, convection in 2017 was preconditioned by the 2016 event, and triggered by an anomalous salt advection within the upper 100 m.

### 3.2.3. Impact of the Vertical Convective Mixing Scheme

To understand the role of the vertical convective mixing scheme adopted in the model in shaping polynya formation and convective activity, we repeat the sensitivity run with 20% precipitation using the different convective mixing schemes (EVD, EDMF and $EDMF_ρ$) after the spin-up period. Since the vertical convective mixing scheme is the only change in the model setup, any differences must relate to the choice of scheme. We focus on the differences in the maximum MLD, the stratification below 50 m beneath the MLD (hereafter $MLD_{+50}$), and relation of $N^2$ to the velocity of the convective plumes in 2016. We use the $MLD_{+50}$ threshold, to exclude mixing within the often highly turbulent mixed layer.

In all runs, a maximum MLD deeper than 800 m is found in 2016 (Figures 8a–8d), but the expansiveness of the area with such deep mixed layers differs between simulations. Using EVD leads to deep MLDs at the northern flank of Maud Rise, as well as toward the south (at 4°W, 68.5°S) in small, confined regions (Figure 8a). The maximum MLD patterns are identical between EDMF and $EDMF_α$ (not shown). Changing the thermal expansion and haline contraction coefficients to non-fixed values within the scheme has no impact on the deep MLD properties. Using $EDMF_ρ$ (Figure 8c) produces the largest region of deep MLDs amongst all runs. This results from changing the potential density to in situ density within the scheme. Including the frictional term (Equation 3), which is related to the lateral entrainment (EDMF$_e$, Figure 8d) leads to a smaller region of deep MLDs. The lateral entrainment term (Equation 4) is rather dominated by the turbulence at the plume-environment interface than the background entrainment $ε_g$. The turbulence term is on average for 2016 around 0.86, while the background entrainment is kept constant with 0.001. Further, it is notable that all EDMF schemes produce a clearer maximum MLD patterns, while with EVD the region south of Maud Rise (which host MLDs shallower than 200 m in all EDMF runs) exhibits MLDs of about 500 m. This difference in convective behavior likely stems from the distinct ways in which the EVD and EDMF schemes represent convective mixing. EDMF works on a sub-grid scale and penetrates through the water column, while EVD mixes neighboring grid cells vertically. In the runs using EVD and 10% and 0% precipitation (Figures 8e and 8f), the maximum MLD reaches below 800 m depth, localized on top of Maud Rise and its flanks. Additionally, deep MLDs (>500 m) can be found north of Astrid Ridge and south of Maud Rise. The reduction from 10% to 0% precipitation leads to an increased area with MLD >800 m and an overall increase in the MLD, but with consistent spatial patterns. To conclude, all tested EDMF schemes are able to produce localized patterns of deep convection that only impact specific regions, while the commonly used EVD generates weaker convection affecting a larger region.

The stratification in the deeper water column is also impacted by the scheme choice. (Monthly means of 2016 for temperature and salinity for the different convective mixing schemes can be found in the Figure S4 in Supporting Information S1.) Therefore, the seasonal-horizontal mean stratification for the Taylor Cap and Flanks (Figure 1 blue and orange encircled area) below $MLD_{+50}$ is estimated for each depth level from the daily output for the four seasons in 2016: summer (JFM), autumn (AMJ), winter (JAS) and spring (OND) (Figure 9). In summer, the distribution of $N^2$ below $MLD_{+50}$ in the Taylor Cap exhibits the largest range of $N^2$ values when using EVD. The EDMF-produced distribution resembles that for EVD, but the maximum (minimum) $N^2$ is smaller (larger). The smallest values of $N^2$ are found when using $EDMF_e$. The distribution of stratification at the Flanks is almost identical between the schemes. The median for both regions is just above 0. In autumn, the largest $N^2$ values for the Taylor Cap are found when using EDMF, and smallest for $EDMF_ρ$. The same patterns are found at the Flanks. The weak stratification at the Flanks and Taylor Cap in $EDMF_ρ$ is likely affected by the onset of deep convection in late June, when the median is above 0. The stratification of the Taylor Cap is affected by the intrusions of colder
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Maximum mixed layer depth in 2016 of the 20% precipitation sensitivity run using (a) enhanced vertical diffusion (EVD) scheme, (b–d) the different modifications of the Eddy-Diffusivity Mass-Flux scheme. (e) 10% precipitation and (f) 0% precipitation using EVD scheme. Gray lines indicate the 2,500, 3,500 and 5,000 m isobaths and the white contour encircles the area of the Taylor Cap.</figcaption></figure>

and fresher surface waters in June 2016 (Figure S5 in Supporting Information S1), which weaken the stratification. In winter, $N^2$ values are one order of magnitude smaller than in any other season. Here, the median for all schemes and both regions is larger than $0.3 \cdot 10^{-6} s^{-2}$. It is notable that the negative tail in the $N^2$ distribution is most pronounced when using EDMF and $EDMF_ρ$. The minimum $N^2$ is smaller when using $EDMF_ε$ compared to EVD. In spring, the largest values of $N^2$ in the Taylor Cap are found with $EDMF_ρ$, while the smallest occur when using EVD. Results from EDMF and $EDMF_ε$ are alike. A similar behavior is found at the Flanks. For both regions, the median is above 0 for all schemes used.

After assessing the variability in stratification below $MLD_{+50}$, we now investigate where in this part of the water column instabilities ($N^2 < 0$) at the Flanks and the Taylor Cap occur, for the runs using EVD and $EDMF_ε$ (Figure 10a). When using EVD, most instabilities occur close to the MLD, and their incidence decreases with depth at both the Flanks and Taylor Cap. When using $EDMF_ε$, a decrease in instability occurrence with depth is also observed, but is enhanced compared to the EVD run. Further, a difference between the Flanks and the Taylor Cap is noted, with fewer instability events occurring in the Taylor Cap than at the Flanks. Finally, we seek to relate the instabilities to the EDMF scheme, and for this purpose we inspect the plume velocity $w_p > 0$ and $N^2 < 0$ for the $EDMF_ε$ in the defined region of the Taylor Cap and Flanks in Figures 10b and 10c. Throughout most of the water column in both analyzed regions more grid cells have a $w_p > 0$ than $N^2 < 0$. An exception is for the depth level
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Maximum mixed layer depth in 2016 of the 20% precipitation sensitivity run using (a) enhanced vertical diffusion (EVD) scheme, (b–d) the different modifications of the Eddy-Diffusivity Mass-Flux scheme. (e) 10% precipitation and (f) 0% precipitation using EVD scheme. Gray lines indicate the 2,500, 3,500 and 5,000 m isobaths and the white contour encircles the area of the Taylor Cap.</figcaption></figure>

and fresher surface waters in June 2016 (Figure S5 in Supporting Information S1), which weaken the stratification. In winter, $N^2$ values are one order of magnitude smaller than in any other season. Here, the median for all schemes and both regions is larger than $0.3 \cdot 10^{-6} s^{-2}$. It is notable that the negative tail in the $N^2$ distribution is most pronounced when using EDMF and $EDMF_ρ$. The minimum $N^2$ is smaller when using $EDMF_ε$ compared to EVD. In spring, the largest values of $N^2$ in the Taylor Cap are found with $EDMF_ρ$, while the smallest occur when using EVD. Results from EDMF and $EDMF_ε$ are alike. A similar behavior is found at the Flanks. For both regions, the median is above 0 for all schemes used.

After assessing the variability in stratification below $MLD_{+50}$, we now investigate where in this part of the water column instabilities ($N^2 < 0$) at the Flanks and the Taylor Cap occur, for the runs using EVD and $EDMF_ε$ (Figure 10a). When using EVD, most instabilities occur close to the MLD, and their incidence decreases with depth at both the Flanks and Taylor Cap. When using $EDMF_ε$, a decrease in instability occurrence with depth is also observed, but is enhanced compared to the EVD run. Further, a difference between the Flanks and the Taylor Cap is noted, with fewer instability events occurring in the Taylor Cap than at the Flanks. Finally, we seek to relate the instabilities to the EDMF scheme, and for this purpose we inspect the plume velocity $w_p > 0$ and $N^2 < 0$ for the $EDMF_ε$ in the defined region of the Taylor Cap and Flanks in Figures 10b and 10c. Throughout most of the water column in both analyzed regions more grid cells have a $w_p > 0$ than $N^2 < 0$. An exception is for the depth level
<figure><img src="image_12.png" /></figure>

**Figure 9.** Distribution of temporal and horizontal mean $N^2$ in the water column below $MLD_{+50}$ for the Maud Rise Flanks (left hand side) and Taylor Cap (right hand side) for the different convective mixing schemes in panel (a) Summer (JFM), (b) Autumn (AMJ), (c) Winter (JAS) and (d) Spring (OND) in 2016. Note the different order of magnitude in panel (c). Dots (stars) indicate the median (mean) of the distribution.

<figure><img src="image_13.png" /></figure>

**Figure 10.** (a) Number of grid cells with a value of $N^2 < 0$ for the Flanks and Taylor Cap using enhanced vertical diffusion (circles) and $EDMF_\epsilon$ (triangles). Relation between grid cells with $w_p > 0$ and $N^2 < 0$ for (b) the Taylor Cap and (c) the Flanks. The dotted line in panel (b, c) indicates equal grid cells for $N^2 < 0$ and $w_p > 0$.
<figure><img src="image_11.png" /></figure>

**Figure 9.** Distribution of temporal and horizontal mean $N^2$ in the water column below $MLD_{+50}$ for the Maud Rise Flanks (left hand side) and Taylor Cap (right hand side) for the different convective mixing schemes in panel (a) Summer (JFM), (b) Autumn (AMJ), (c) Winter (JAS) and (d) Spring (OND) in 2016. Note the different order of magnitude in panel (c). Dots (stars) indicate the median (mean) of the distribution.

<figure><img src="image_12.png" /></figure>

**Figure 10.** (a) Number of grid cells with a value of $N^2 < 0$ for the Flanks and Taylor Cap using enhanced vertical diffusion (circles) and $EDMF_\epsilon$ (triangles). Relation between grid cells with $w_p > 0$ and $N^2 < 0$ for (b) the Taylor Cap and (c) the Flanks. The dotted line in panel (b, c) indicates equal grid cells for $N^2 < 0$ and $w_p > 0$.
close to the bottom, where only active $N^2 < 0$ can be found due to bottom mixing. This difference in the upper ocean is mainly due to the scheme's ability to connect several vertical grid cells. This allows us to conclude that $N^2 < 0$ can be used as an indicator of the EDMF scheme's production of a convective plume. In conclusion, EDMF generates more mixing at greater depth than the EVD.

## 4. Discussion and Conclusion

In this study, we were able to reproduce the MRPs observed in 2016 and 2017 in a regional ocean-ice configuration with a forced atmosphere, and we highlighted the importance of vertical convective mixing schemes and surface freshwater forcing in determining the extent and nature of polynya formation and associated deep convection. In a series of numerical experiments, we tested different vertical convective mixing schemes, reduced freshwater forcings, and a combination to study their impact on the production and characteristics of MRPs.

Simulations using EVD, EDMF and $EDMF_c$ for the period 2007–2017 with unperturbed precipitation did not lead to a polynya formation, but were used as non-polynya runs and compared to observational climatologies (Figure 2). From the differences between observations and simulations, it was shown that the convective mixing scheme impacts the MLD and the ocean properties below it, and that the modifications applied to the EDMF scheme led to a better representation of MLD and smaller biases in the upper 600 m. Overall, the simulations showed a fresh and cold bias, especially at the surface. The bias in the simulations resulted from excessive ice transport from the open boundary conditions into the domain (Table 2), which led to an overly strong stratification that inhibited convection. The comparison to observations is useful but can contain large uncertainties related to sampling biases due to the sparse availability of observations in the region.

Based on the non-polynya simulations, a set of sensitivity runs was conducted to reduce the freshening by sea-ice transport from the open boundaries. Accordingly, the precipitation was reduced compared to the original precipitation while using $EDMF_c$ (Figure 3). The first polynya signatures were detected in runs using only 30% of the original precipitation. In runs with 10% or less, polynyas occurred in more years than observed. In our simulations, the reduction of precipitation compensated for the unrealistically high convergence of sea-ice in the model domain. The sensitivity runs show that simulations are very sensitive to the freshwater forcing, and the polynya activity is easily over- or underestimated. That freshwater forcing plays a key role with respect to convection events has been shown before, but most often the freshwater forcing was increased to suppress convection events and spurious polynya openings (Kjellsson et al., 2015; Stoessel et al., 2015). Additionally to the ice compensation, the lack of observations of precipitation in the Southern Ocean impacts the accuracy of all precipitation products (Siems et al., 2022). This can lead to a substantial over- or underestimation of precipitation in JRA55-do version 1.5, which is used as surface forcing in our simulations.

For the analysis of the polynya formation processes, we investigated the properties of the 20% precipitation sensitivity run. This run showed a high coherence with observations regarding polynya location, even though a time difference between model and observations exists. Our analysis showed that, in 2016, the polynya resulted from an instability between 200 and 250 m in combination with a reduced ambient stratification (Figure 7). This is supported by the works of Akitomo (2006) and McPhee (2000), who showed that the flanks of Maud Rise are susceptible to thermobaric instabilities. Similar to Akitomo (2006), we found that the water column at the flanks is more unstable than on top of Maud Rise (Figure 10a). In 2017, a positive salt anomaly preceding the onset of deep convection was detected in the upper 100 m, which led to an overturning of the water column. This salt anomaly in the surface has been found in other model and observational studies of the 2017 Maud Rise Polynya (Campbell et al., 2019; Kurtakoti et al., 2018). Our analysis highlighted that the deep convection of 2016 preconditioned the water column for the opening in 2017. This enabled the larger and more persistent polynya of 2017, in agreement with the findings of Campbell et al. (2019).

In our last experiment, we repeated the 20% precipitation run with the different vertical convective mixing schemes. In general, the use of EDMF led to more defined patterns of deep mixed layers, and suppressed extensive convection throughout the entire domain (Figure 8). Further, EDMF enhanced mixing at depth (Figure 10). Enhancing vertical mixing was pointed out in previous studies as a possible way to better represent the MLD and deep convection (Heuzé et al., 2015; Kjellsson et al., 2015). Additionally, Kjellsson et al. (2015) and Stoessel et al. (2015) highlighted the importance of freshwater forcing in the control of convection in the Southern Ocean. Our approach of combining freshwater forcing and vertical mixing points to the importance of finding an optimal balance between freshwater forcing and vertical mixing. This proves to be a challenge, due to the
close to the bottom, where only active $N^2 < 0$ can be found due to bottom mixing. This difference in the upper ocean is mainly due to the scheme's ability to connect several vertical grid cells. This allows us to conclude that $N^2 < 0$ can be used as an indicator of the EDMF scheme's production of a convective plume. In conclusion, EDMF generates more mixing at greater depth than the EVD.

## 4. Discussion and Conclusion

In this study, we were able to reproduce the MRPs observed in 2016 and 2017 in a regional ocean-ice configuration with a forced atmosphere, and we highlighted the importance of vertical convective mixing schemes and surface freshwater forcing in determining the extent and nature of polynya formation and associated deep convection. In a series of numerical experiments, we tested different vertical convective mixing schemes, reduced freshwater forcings, and a combination to study their impact on the production and characteristics of MRPs.

Simulations using EVD, EDMF and $EDMF_c$ for the period 2007–2017 with unperturbed precipitation did not lead to a polynya formation, but were used as non-polynya runs and compared to observational climatologies (Figure 2). From the differences between observations and simulations, it was shown that the convective mixing scheme impacts the MLD and the ocean properties below it, and that the modifications applied to the EDMF scheme led to a better representation of MLD and smaller biases in the upper 600 m. Overall, the simulations showed a fresh and cold bias, especially at the surface. The bias in the simulations resulted from excessive ice transport from the open boundary conditions into the domain (Table 2), which led to an overly strong stratification that inhibited convection. The comparison to observations is useful but can contain large uncertainties related to sampling biases due to the sparse availability of observations in the region.

Based on the non-polynya simulations, a set of sensitivity runs was conducted to reduce the freshening by sea-ice transport from the open boundaries. Accordingly, the precipitation was reduced compared to the original precipitation while using $EDMF_c$ (Figure 3). The first polynya signatures were detected in runs using only 30% of the original precipitation. In runs with 10% or less, polynyas occurred in more years than observed. In our simulations, the reduction of precipitation compensated for the unrealistically high convergence of sea-ice in the model domain. The sensitivity runs show that simulations are very sensitive to the freshwater forcing, and the polynya activity is easily over- or underestimated. That freshwater forcing plays a key role with respect to convection events has been shown before, but most often the freshwater forcing was increased to suppress convection events and spurious polynya openings (Kjellsson et al., 2015; Stoessel et al., 2015). Additionally to the ice compensation, the lack of observations of precipitation in the Southern Ocean impacts the accuracy of all precipitation products (Siems et al., 2022). This can lead to a substantial over- or underestimation of precipitation in JRA55-do version 1.5, which is used as surface forcing in our simulations.

For the analysis of the polynya formation processes, we investigated the properties of the 20% precipitation sensitivity run. This run showed a high coherence with observations regarding polynya location, even though a time difference between model and observations exists. Our analysis showed that, in 2016, the polynya resulted from an instability between 200 and 250 m in combination with a reduced ambient stratification (Figure 7). This is supported by the works of Akitomo (2006) and McPhee (2000), who showed that the flanks of Maud Rise are susceptible to thermobaric instabilities. Similar to Akitomo (2006), we found that the water column at the flanks is more unstable than on top of Maud Rise (Figure 10a). In 2017, a positive salt anomaly preceding the onset of deep convection was detected in the upper 100 m, which led to an overturning of the water column. This salt anomaly in the surface has been found in other model and observational studies of the 2017 Maud Rise Polynya (Campbell et al., 2019; Kurtakoti et al., 2018). Our analysis highlighted that the deep convection of 2016 preconditioned the water column for the opening in 2017. This enabled the larger and more persistent polynya of 2017, in agreement with the findings of Campbell et al. (2019).

In our last experiment, we repeated the 20% precipitation run with the different vertical convective mixing schemes. In general, the use of EDMF led to more defined patterns of deep mixed layers, and suppressed extensive convection throughout the entire domain (Figure 8). Further, EDMF enhanced mixing at depth (Figure 10). Enhancing vertical mixing was pointed out in previous studies as a possible way to better represent the MLD and deep convection (Heuzé et al., 2015; Kjellsson et al., 2015). Additionally, Kjellsson et al. (2015) and Stoessel et al. (2015) highlighted the importance of freshwater forcing in the control of convection in the Southern Ocean. Our approach of combining freshwater forcing and vertical mixing points to the importance of finding an optimal balance between freshwater forcing and vertical mixing. This proves to be a challenge, due to the
complexity and parameter options related to vertical mixing, as well as to the interplay with other factors, for example, model resolution.

Here, for the first time, an ocean-ice model forced by an atmospheric reanalysis was able to reproduce the MRP events of 2016 and 2017. The 2016 MRP was triggered by a thermobaric instability below the MLD, between 200 and 250 m depth in combination with a reduced ambient stratification. The aforementioned depth range encompasses the interface between WDW and the cold surface layer. The thermobaric instability was enabled in the modified vertical convective mixing scheme $EDMF_c$. The 2017 polynya resulted from a preconditioned weak stratification combined with a surface salinity advection event. In both cases, the EDMF schemes led to localized deep convection. Based on our results, we recommend generalizing the use of penetrative convection schemes such as EDMF in all OGCMs, since it helps to reproduce more realistic convective plumes and events, and it also improves the representation of the MLD in non-convective periods.

## Data Availability Statement

The EN4.2.2 data set (Met Office Hadley Centre, 2022) is available from Met Office Hadley Center at https://www.metoffice.gov.uk/hadobs/en4/. The CTD data from MaudNESS (McPhee, 2005) can be found at: https://doi.org/10.7284/152517 and the SO-CHIC CTD data (Steiger et al., 2022) can be found at: https://doi.org/10.17882/95314. The SIC satellite data from AMSR-E (Sreen et al., 2008) and AMSR2 (Melsheimer & Sreen, 2019) are available from the University of Bremen: https://seaice.uni-bremen.de/data. ERA5 reanalysis data (Hersbach et al., 2023) is available at the European Center for Medium Range Weather Forecast (ECMWF) via: https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels. For the inlay in Figure 1 the bathymetry data from GEBCO (GEBCO Compilation Group, 2023) is used. All computer code used to generate the model output and computer code to reproduce the plotted figures are available in an open-access repository: https://github.com/bguelk/JAMES_2024_MRP_modelling (Gülk, 2024).

## References

* Akitomo, K. (1999a). Open-ocean deep convection due to thermobaricity: 1. Scaling argument. *Journal of Geophysical Research*, 104(C3), 5225-5234. https://doi.org/10.1029/1998jc900058

* Akitomo, K. (1999b). Open-ocean deep convection due to thermobaricity: 2. Numerical experiments. *Journal of Geophysical Research*, 104(C3), 5235-5249. https://doi.org/10.1029/1998jc900062

* Akitomo, K. (2006). Thermobaric deep convection, baroclinic instability, and their roles in vertical heat transport around Maud Rise in the Weddell Sea. *Journal of Geophysical Research*, 111(C9). https://doi.org/10.1029/2005jc003284

* Akitomo, K., Awaji, T., & Imasato, N. (1995). Open-ocean deep convection in the Weddell Sea: Two-dimensional numerical experiments with a nonhydrostatic model. *Deep Sea Research Part I: Oceanographic Research Papers*, 42(1), 53-73. https://doi.org/10.1016/0967-0637(94)00035-q

* Bernardello, R., Marinov, I., Palter, J. B., Galbraith, E. D., & Sarmiento, J. L. (2014). Impact of Weddell Sea deep convection on natural and anthropogenic carbon in a climate model. *Geophysical Research Letters*, 41(20), 7262-7269. https://doi.org/10.1002/2014gl061313

* Campbell, E. C., Wilson, E. A., Moore, G. K., Riser, S. C., Brayton, C. E., Mazloff, M. R., & Talley, L. D. (2019). Antarctic offshore polynyas linked to southern hemisphere climate anomalies. *Nature*, 570(7761), 319-325. https://doi.org/10.1038/s41586-019-1294-0

* Carsey, F. (1980). Microwave observation of the Weddell polynya. *Monthly Weather Review*, 108(12), 2032-2044. https://doi.org/10.1175/1520-0493(1980)108<2032:mootwp>2.0.co;2

* Cheon, W. G., & Gordon, A. L. (2019). Open-ocean polynyas and deep convection in the southern ocean. *Scientific Reports*, 9(1), 1-9. https://doi.org/10.1038/s41598-019-43466-2

* Cheon, W. G., Lee, S.-K., Gordon, A. L., Liu, Y., Cho, C.-B., & Park, J. J. (2015). Replicating the 1970s' Weddell polynya using a coupled ocean-sea ice model with reanalysis surface flux fields. *Geophysical Research Letters*, 42(13), 5411-5418. https://doi.org/10.1002/2015gl064364

* Cisewski, B., Strass, V. H., & Leach, H. (2011). Circulation and transport of water masses in the Lazarev Sea, Antarctica, during summer and winter 2006. *Deep Sea Research Part I: Oceanographic Research Papers*, 58(2), 186-199. https://doi.org/10.1016/j.dsr.2010.12.001

* Comiso, J., & Gordon, A. (1987). Recurring polynyas over the cosmonaut sea and the Maud Rise. *Journal of Geophysical Research*, 92(C3), 2819-2833. https://doi.org/10.1029/jc092ic03p02819

* De Steur, L., Holland, D., Muench, R., & McPhee, M. G. (2007). The warm-water "halo" around Maud Rise: Properties, dynamics and impact. *Deep Sea Research Part I: Oceanographic Research Papers*, 54(6), 871-896. https://doi.org/10.1016/j.dsr.2007.03.009

* Dufour, C. O., Morrison, A. K., Griffies, S. M., Frenger, I., Zanowski, H., & Winton, M. (2017). Preconditioning of the Weddell Sea polynya by the ocean mesoscale and dense water overflows. *Journal of Climate*, 30(19), 7719-7737. https://doi.org/10.1175/jcli-d-16-0586.1

* Francis, D., Eayrs, C., Cuesta, J., & Holland, D. (2019). Polar cyclones at the origin of the reoccurrence of the Maud Rise Polynya in austral winter 2017. *Journal of Geophysical Research: Atmospheres*, 124(10), 5251-5267. https://doi.org/10.1029/2019jd030618

* Francis, D., Mattingly, K. S., Temimi, M., Massom, R., & Heil, P. (2020). On the crucial role of atmospheric rivers in the two major Weddell polynya events in 1973 and 2017 in Antarctica. *Science Advances*, 6(46), eabc2695. https://doi.org/10.1126/sciadv.abc2695

* GEBCO Compilation Group. (2023). Gebco 2023 grid [Dataset]. https://doi.org/10.5285/f98b053b-0cbc-6c23-e053-6c86abc0af7b

* Giordani, H., Bourdallé-Badie, R., & Madec, G. (2020). An eddy-diffusivity mass-flux parameterization for modeling oceanic convection. *Journal of Advances in Modeling Earth Systems*, 12(9), e2020MS002078. https://doi.org/10.1029/2020ms002078

* Good, S. A., Martin, M. J., & Rayner, N. A. (2013). En4: Quality controlled ocean temperature and salinity profiles and monthly objective analyses with uncertainty estimates. *Journal of Geophysical Research: Oceans*, 118(12), 6704-6716. https://doi.org/10.1002/2013jc009067

## Acknowledgments

The authors thank the three anonymous reviewers and the editor for their useful comments and suggestions. BG, FR, ACNG and GM acknowledge support from European Union's Horizon 2020 research and innovation programme under Grant 821001 (SO-CHIC). The computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) and the Swedish National Infrastructure for Computing (SNIC) at NSC partially funded by the Swedish Research Council through Grant 2022-06725 and 2018-05973.
complexity and parameter options related to vertical mixing, as well as to the interplay with other factors, for example, model resolution.

Here, for the first time, an ocean-ice model forced by an atmospheric reanalysis was able to reproduce the MRP events of 2016 and 2017. The 2016 MRP was triggered by a thermobaric instability below the MLD, between 200 and 250 m depth in combination with a reduced ambient stratification. The aforementioned depth range encompasses the interface between WDW and the cold surface layer. The thermobaric instability was enabled in the modified vertical convective mixing scheme $EDMF_c$. The 2017 polynya resulted from a preconditioned weak stratification combined with a surface salinity advection event. In both cases, the EDMF schemes led to localized deep convection. Based on our results, we recommend generalizing the use of penetrative convection schemes such as EDMF in all OGCMs, since it helps to reproduce more realistic convective plumes and events, and it also improves the representation of the MLD in non-convective periods.

## Data Availability Statement

The EN4.2.2 data set (Met Office Hadley Centre, 2022) is available from Met Office Hadley Center at https://www.metoffice.gov.uk/hadobs/en4/. The CTD data from MaudNESS (McPhee, 2005) can be found at: https://doi.org/10.7284/152517 and the SO-CHIC CTD data (Steiger et al., 2022) can be found at: https://doi.org/10.17882/95314. The SIC satellite data from AMSR-E (Sreen et al., 2008) and AMSR2 (Melsheimer & Sreen, 2019) are available from the University of Bremen: https://seaice.uni-bremen.de/data. ERA5 reanalysis data (Hersbach et al., 2023) is available at the European Center for Medium Range Weather Forecast (ECMWF) via: https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels. For the inlay in Figure 1 the bathymetry data from GEBCO (GEBCO Compilation Group, 2023) is used. All computer code used to generate the model output and computer code to reproduce the plotted figures are available in an open-access repository: https://github.com/bguelk/JAMES_2024_MRP_modelling (Gülk, 2024).

## References

* Akitomo, K. (1999a). Open-ocean deep convection due to thermobaricity: 1. Scaling argument. *Journal of Geophysical Research*, 104(C3), 5225–5234. https://doi.org/10.1029/1998jc900058

* Akitomo, K. (1999b). Open-ocean deep convection due to thermobaricity: 2. Numerical experiments. *Journal of Geophysical Research*, 104(C3), 5235–5249. https://doi.org/10.1029/1998jc900062

* Akitomo, K. (2006). Thermobaric deep convection, baroclinic instability, and their roles in vertical heat transport around Maud Rise in the Weddell Sea. *Journal of Geophysical Research*, 111(C9). https://doi.org/10.1029/2005jc003284

* Akitomo, K., Awaji, T., & Imasato, N. (1995). Open-ocean deep convection in the Weddell Sea: Two-dimensional numerical experiments with a nonhydrostatic model. *Deep Sea Research Part I: Oceanographic Research Papers*, 42(1), 53–73. https://doi.org/10.1016/0967-0637(94)00035-q

* Bernardello, R., Marinov, I., Palter, J. B., Galbraith, E. D., & Sarmiento, J. L. (2014). Impact of Weddell Sea deep convection on natural and anthropogenic carbon in a climate model. *Geophysical Research Letters*, 41(20), 7262–7269. https://doi.org/10.1002/2014gl061313

* Campbell, E. C., Wilson, E. A., Moore, G. K., Riser, S. C., Brayton, C. E., Mazloff, M. R., & Talley, L. D. (2019). Antarctic offshore polynyas linked to southern hemisphere climate anomalies. *Nature*, 570(7761), 319–325. https://doi.org/10.1038/s41586-019-1294-0

* Carsey, F. (1980). Microwave observation of the Weddell polynya. *Monthly Weather Review*, 108(12), 2032–2044. https://doi.org/10.1175/1520-0493(1980)108<2032:mootwp>2.0.co;2

* Cheon, W. G., & Gordon, A. L. (2019). Open-ocean polynyas and deep convection in the southern ocean. *Scientific Reports*, 9(1), 1–9. https://doi.org/10.1038/s41598-019-43466-2

* Cheon, W. G., Lee, S.-K., Gordon, A. L., Liu, Y., Cho, C.-B., & Park, J. J. (2015). Replicating the 1970s' Weddell polynya using a coupled ocean-sea ice model with reanalysis surface flux fields. *Geophysical Research Letters*, 42(13), 5411–5418. https://doi.org/10.1002/2015gl064364

* Cisewski, B., Strass, V. H., & Leach, H. (2011). Circulation and transport of water masses in the Lazarev Sea, Antarctica, during summer and winter 2006. *Deep Sea Research Part I: Oceanographic Research Papers*, 58(2), 186–199. https://doi.org/10.1016/j.dsr.2010.12.001

* Comiso, J., & Gordon, A. (1987). Recurring polynyas over the cosmonaut sea and the Maud Rise. *Journal of Geophysical Research*, 92(C3), 2819–2833. https://doi.org/10.1029/jc092ic03p02819

* De Steur, L., Holland, D., Muench, R., & McPhee, M. G. (2007). The warm-water “halo” around Maud Rise: Properties, dynamics and impact. *Deep Sea Research Part I: Oceanographic Research Papers*, 54(6), 871–896. https://doi.org/10.1016/j.dsr.2007.03.009

* Dufour, C. O., Morrison, A. K., Griffies, S. M., Frenger, I., Zanowski, H., & Winton, M. (2017). Preconditioning of the Weddell Sea polynya by the ocean mesoscale and dense water overflows. *Journal of Climate*, 30(19), 7719–7737. https://doi.org/10.1175/jcli-d-16-0586.1

* Francis, D., Eayrs, C., Cuesta, J., & Holland, D. (2019). Polar cyclones at the origin of the reoccurrence of the Maud Rise Polynya in austral winter 2017. *Journal of Geophysical Research: Atmospheres*, 124(10), 5251–5267. https://doi.org/10.1029/2019jd030618

* Francis, D., Mattingly, K. S., Temimi, M., Massom, R., & Heil, P. (2020). On the crucial role of atmospheric rivers in the two major Weddell polynya events in 1973 and 2017 in Antarctica. *Science Advances*, 6(46), eabc2695. https://doi.org/10.1126/sciadv.abc2695

* GEBCO Compilation Group. (2023). Gebco 2023 grid [Dataset]. https://doi.org/10.5285/f98b053b-0cbc-6c23-e053-6c86abc0af7b

* Giordani, H., Bourdallé-Badie, R., & Madec, G. (2020). An eddy-diffusivity mass-flux parameterization for modeling oceanic convection. *Journal of Advances in Modeling Earth Systems*, 12(9), e2020MS002078. https://doi.org/10.1029/2020ms002078

* Good, S. A., Martin, M. J., & Rayner, N. A. (2013). En4: Quality controlled ocean temperature and salinity profiles and monthly objective analyses with uncertainty estimates. *Journal of Geophysical Research: Oceans*, 118(12), 6704–6716. https://doi.org/10.1002/2013jc009067

## Acknowledgments

The authors thank the three anonymous reviewers and the editor for their useful comments and suggestions. BG, FR, ACNG and GM acknowledge support from European Union's Horizon 2020 research and innovation programme under Grant 821001 (SO-CHIC). The computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) and the Swedish National Infrastructure for Computing (SNIC) at NSC partially funded by the Swedish Research Council through Grant 2022-06725 and 2018-05973.
Gordon, A. L. (1982). Weddell deep water variability. *Journal of Marine Research*, 40, 199-217.

Gordon, A. L., & Huber, B. A. (1990). Southern ocean winter mixed layer. *Journal of Geophysical Research*, 95(C7), 11655-11672. https://doi.org/10.1029/jc095ic07p11655

Gregory, D. (2001). Estimation of entrainment rate in simple models of convective clouds. *Quarterly Journal of the Royal Meteorological Society*, 127(571), 53-72. https://doi.org/10.1256/smsqi.57103

Gülk, B. (2024). EDMF routines, Maud12 parameters and data analysis [Software]. https://github.com/bguelk/JAMES_2024_MRP_modelling

Gülk, B., Roquet, F., Naveira Garabato, A. C., Narayanan, A., Rousset, C., & Madec, G. (2023). Variability and remote controls of the warm-water halo and Taylor Cap at Maud Rise. *Journal of Geophysical Research: Oceans*, 128(7), e2022JC019517. https://doi.org/10.1029/2022jc019517

Hersbach, H., Bell, B., Berrisford, P., Biavati, G., Horányi, A., Muñoz Sabater, J., et al. (2023). Era5 hourly data on single levels from 1940 to present [Dataset]. https://doi.org/10.24381/cds.adbb2d47

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., et al. (2020). The ERA5 global reanalysis. *Quarterly Journal of the Royal Meteorological Society*, 146(730), 1999-2049. https://doi.org/10.1002/qj.3803

Heuzé, C., Ridley, J., Calvert, D., Stevens, D., & Heywood, K. (2015). Increasing vertical mixing to reduce southern ocean deep convection in NEMO3. 4. *Geoscientific Model Development*, 8(10), 3119-3130. https://doi.org/10.5194/gmd-8-3119-2015

Holland, D. (2001). Explaining the Weddell polynya-a large ocean eddy shed at Maud Rise. *Science*, 292(5522), 1697-1700. https://doi.org/10.1126/science.1059322

Kjellsson, J., Holland, P. R., Marshall, G. J., Mathiot, P., Aksenov, Y., Coward, A. C., et al. (2015). Model sensitivity of the Weddell and Ross seas, Antarctica, to vertical mixing and freshwater forcing. *Ocean Modelling*, 94, 141-152. https://doi.org/10.1016/j.ocemod.2015.08.003

Kurtakoti, P., Veneziani, M., Stössel, A., & Weijer, W. (2018). Preconditioning and formation of Maud Rise Polynyas in a high-resolution earth system model. *Journal of Climate*, 31(23), 9659-9678. https://doi.org/10.1175/jcli-d-18-0392.1

Kurtakoti, P., Veneziani, M., Stössel, A., Weijer, W., & Maltrud, M. (2021). On the generation of Weddell Sea polynyas in a high-resolution earth system model. *Journal of Climate*, 34(7), 2491-2510. https://doi.org/10.1175/jcli-d-20-0229.1

Leach, H., Strass, V., & Cisewski, B. (2011). Modification by lateral mixing of the warm deep water entering the Weddell Sea in the Maud Rise region. *Ocean Dynamics*, 61(1), 51-68. https://doi.org/10.1007/s10236-010-0342-y

Lellouche, J.-M., Greiner, E., Bourdallé-Badie, R., Garric, G., Melet, A., Drévillon, M., et al. (2021). The Copernicus global 1/12° oceanic and sea ice GLORYS12 reanalysis. *Frontiers in Earth Science*, 9. https://doi.org/10.3389/feart.2021.698876

Lindsay, R., Holland, D., & Woodgate, R. (2004). Halo of low ice concentration observed over the Maud Rise seamount. *Geophysical Research Letters*, 31(13). https://doi.org/10.1029/2004gl019831

Madec, G., & NEMO System Team. (2019). NEMO ocean engine. *Scientific notes of climate modelling center*. Institut Pierre-Simon Laplace (IPSL). https://doi.org/10.5281/ZENODO.1464816

Marshall, J., & Schott, F. (1999). Open-ocean convection: Observations, theory, and models. *Reviews of Geophysics*, 37(1), 1-64. https://doi.org/10.1029/98rg02739

Martinson, D. G., & Iannuzzi, R. A. (1998). Antarctic ocean-ice interaction: Implications from ocean bulk property distributions in the Weddell Gyre. *Antarctic sea ice: physical processes, interactions and variability*, 74, 243-271. https://doi.org/10.1029/ar074p0243

Mazloff, M. R., Heimbach, P., & Wunsch, C. (2010). An eddy-permitting southern ocean state estimate. *Journal of Physical Oceanography*, 40(5), 880-899. https://doi.org/10.1175/2009jpo4236.1

McPhee, M. G. (2000). Marginal thermobaric stability in the ice-covered upper ocean over Maud Rise. *Journal of Physical Oceanography*, 30(11), 2710-2722. https://doi.org/10.1175/1520-0485(2000)030<2710:mtsiti>2.0.co;2

McPhee, M. G. (2003). Is thermobaricity a major factor in southern ocean ventilation? *Antarctic Science*, 15(1), 153-160. https://doi.org/10.1017/s0954102003001159

McPhee, M. G. (2005). Ctd (seabird SBE-911+) data as collected during the cruise NBP0506 [Dataset]. *Rolling Deck to Repository (R2R).* https://doi.org/10.7284/152517

Melsheimer, C., & Sreen, G. (2019). AMSR2 ASI sea ice concentration data, Antarctic, version 5.4 (NetCDF) [Dataset]. PANGAEA. https://doi.org/10.1594/PANGAEA.898400

Met Office Hadley Centre. (2022). EN4 quality controlled ocean data [Dataset]. https://www.metoffice.gov.uk/hadobs/en4/

Mohrmann, M., Heuzé, C., & Swart, S. (2021). Southern ocean polynyas in CMIP6 models. *The Cryosphere*, 15(9), 4281-4313. https://doi.org/10.5194/tc-15-4281-2021

NEMO Sea Ice Working Group. (2018). Sea ice modelling integrated initiative (SI3) - The NEMO sea ice engine (No. 31). *Zenodo*. https://doi.org/10.5281/zenodo.1471689

Rheinländer, J. W., Smedsrud, L. H., & Nisancioglu, K. H. (2021). Internal ocean dynamics control the long-term evolution of Weddell Sea polynya activity. *Frontiers in Climate*, 3, 718016. https://doi.org/10.3389/fclim.2021.718016

Roquet, F., Madec, G., McDougall, T. J., & Barker, P. M. (2015). Accurate polynomial expressions for the density and specific volume of seawater using the TEOS-10 standard. *Ocean Modelling*, 90, 29-43. https://doi.org/10.1016/j.ocemod.2015.04.002

Schröder, M., & Fahrbach, E. (1999). On the structure and the transport of the eastern Weddell gyre. *Deep Sea Research Part II: Topical Studies in Oceanography*, 46(1-2), 501-527. https://doi.org/10.1016/s0967-0645(98)00112-x

Siebesma, A. P., Bretherton, C. S., Brown, A., Chlond, A., Cuxart, J., Duynkterke, P. G., et al. (2003). A large eddy simulation intercomparison study of shallow cumulus convection. *Journal of the Atmospheric Sciences*, 60(10), 1201-1219. https://doi.org/10.1175/1520-0469(2003)60<1201:alesis>2.0.co;2

Siems, S. T., Huang, Y., & Manton, M. J. (2022). Southern ocean precipitation: Toward a process-level understanding. *Wiley Interdisciplinary Reviews: Climate Change*, 13(6), e800. https://doi.org/10.1002/wcc.800

Sreen, G., Kaleschke, L., & Heygster, G. (2008). Sea ice remote sensing using AMSR-e 89-GHZ channels. *Journal of Geophysical Research*, 113(C2). https://doi.org/10.1029/2005jc003384

Steiger, N., Sallée, J., Ward, B., Azevedo, A., Binase, Z., Ten Doeschat, A., et al. (2022). CTD observation - So-chic cruise 2022 [Dataset]. *SEANOE*. https://doi.org/10.17882/95314

Stoessel, A., Notz, D., Haumann, F. A., Haak, H., Jungclaus, J., & Mikolajewicz, U. (2015). Controlling high-latitude southern ocean convection in climate models. *Ocean Modelling*, 86, 58-75. https://doi.org/10.1016/j.ocemod.2014.11.008

Tsujino, H., Urakawa, S., Nakano, H., Small, R. J., Kim, W. M., Yeager, S. G., et al. (2018). JRA-55 based surface dataset for driving ocean-sea-ice models (JRA55-DO). *Ocean Modelling*, 130, 79-139. https://doi.org/10.1016/j.ocemod.2018.07.002

Wilson, E. A., Riser, S. C., Campbell, E. C., & Wong, A. P. (2019). Winter upper-ocean stability and ice-ocean feedbacks in the sea ice-covered southern ocean. *Journal of Physical Oceanography*, 49(4), 1099-1117. https://doi.org/10.1175/jpo-d-18-0184.1

Zhou, L., Heuzé, C., & Mohrmann, M. (2022). Early winter triggering of the Maud Rise Polynya. *Geophysical Research Letters*, 49(2), e2021GL096246. https://doi.org/10.1029/2021GL096246
Gordon, A. L. (1982). Weddell deep water variability. *Journal of Marine Research*, 40, 199-217.

Gordon, A. L., & Huber, B. A. (1990). Southern ocean winter mixed layer. *Journal of Geophysical Research*, 95(C7), 11655-11672. https://doi.org/10.1029/jc095ic07p11655

Gregory, D. (2001). Estimation of entrainment rate in simple models of convective clouds. *Quarterly Journal of the Royal Meteorological Society*, 127(571), 53-72. https://doi.org/10.1256/smsqi.57103

Gülk, B. (2024). EDMF routines, Maud12 parameters and data analysis [Software]. https://github.com/bguelk/JAMES_2024_MRP_modelling

Gülk, B., Roquet, F., Naveira Garabato, A. C., Narayanan, A., Rousset, C., & Madec, G. (2023). Variability and remote controls of the warm-water halo and Taylor Cap at Maud Rise. *Journal of Geophysical Research: Oceans*, 128(7), e2022JC019517. https://doi.org/10.1029/2022jc019517

Hersbach, H., Bell, B., Berrisford, P., Biavati, G., Horányi, A., Muñoz Sabater, J., et al. (2023). Era5 hourly data on single levels from 1940 to present [Dataset]. https://doi.org/10.24381/cds.adbb2d47

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., et al. (2020). The ERA5 global reanalysis. *Quarterly Journal of the Royal Meteorological Society*, 146(730), 1999-2049. https://doi.org/10.1002/qj.3803

Heuzé, C., Ridley, J., Calvert, D., Stevens, D., & Heywood, K. (2015). Increasing vertical mixing to reduce southern ocean deep convection in NEMO3. 4. *Geoscientific Model Development*, 8(10), 3119-3130. https://doi.org/10.5194/gmd-8-3119-2015

Holland, D. (2001). Explaining the Weddell polynya-a large ocean eddy shed at Maud Rise. *Science*, 292(5522), 1697-1700. https://doi.org/10.1126/science.1059322

Kjellsson, J., Holland, P. R., Marshall, G. J., Mathiot, P., Aksenov, Y., Coward, A. C., et al. (2015). Model sensitivity of the Weddell and Ross seas, Antarctica, to vertical mixing and freshwater forcing. *Ocean Modelling*, 94, 141-152. https://doi.org/10.1016/j.ocemod.2015.08.003

Kurtakoti, P., Veneziani, M., Stössel, A., & Weijer, W. (2018). Preconditioning and formation of Maud Rise Polynyas in a high-resolution earth system model. *Journal of Climate*, 31(23), 9659-9678. https://doi.org/10.1175/jcli-d-18-0392.1

Kurtakoti, P., Veneziani, M., Stössel, A., Weijer, W., & Maltrud, M. (2021). On the generation of Weddell Sea polynyas in a high-resolution earth system model. *Journal of Climate*, 34(7), 2491-2510. https://doi.org/10.1175/jcli-d-20-0229.1

Leach, H., Strass, V., & Cisewski, B. (2011). Modification by lateral mixing of the warm deep water entering the Weddell Sea in the Maud Rise region. *Ocean Dynamics*, 61(1), 51-68. https://doi.org/10.1007/s10236-010-0342-y

Lellouche, J.-M., Greiner, E., Bourdallé-Badie, R., Garric, G., Melet, A., Drévillon, M., et al. (2021). The Copernicus global 1/12° oceanic and sea ice GLORYS12 reanalysis. *Frontiers in Earth Science*, 9. https://doi.org/10.3389/feart.2021.698876

Lindsay, R., Holland, D., & Woodgate, R. (2004). Halo of low ice concentration observed over the Maud Rise seamount. *Geophysical Research Letters*, 31(13). https://doi.org/10.1029/2004gl019831

Madec, G., & NEMO System Team. (2019). NEMO ocean engine. *Scientific notes of climate modelling center*. Institut Pierre-Simon Laplace (IPSL). https://doi.org/10.5281/ZENODO.1464816

Marshall, J., & Schott, F. (1999). Open-ocean convection: Observations, theory, and models. *Reviews of Geophysics*, 37(1), 1-64. https://doi.org/10.1029/98rg02739

Martinson, D. G., & Iannuzzi, R. A. (1998). Antarctic ocean-ice interaction: Implications from ocean bulk property distributions in the Weddell Gyre. *Antarctic sea ice: physical processes, interactions and variability*, 74, 243-271. https://doi.org/10.1029/ar074p0243

Mazloff, M. R., Heimbach, P., & Wunsch, C. (2010). An eddy-permitting southern ocean state estimate. *Journal of Physical Oceanography*, 40(5), 880-899. https://doi.org/10.1175/2009jpo4236.1

McPhee, M. G. (2000). Marginal thermobaric stability in the ice-covered upper ocean over Maud Rise. *Journal of Physical Oceanography*, 30(11), 2710-2722. https://doi.org/10.1175/1520-0485(2000)030<2710:mtsiti>2.0.co;2

McPhee, M. G. (2003). Is thermobaricity a major factor in southern ocean ventilation? *Antarctic Science*, 15(1), 153-160. https://doi.org/10.1017/s0954102003001159

McPhee, M. G. (2005). Ctd (seabird SBE-911+) data as collected during the cruise NBP0506 [Dataset]. *Rolling Deck to Repository (R2R).* https://doi.org/10.7284/152517

Melsheimer, C., & Sreen, G. (2019). AMSR2 ASI sea ice concentration data, Antarctic, version 5.4 (NetCDF) [Dataset]. PANGAEA. https://doi.org/10.1594/PANGAEA.898400

Met Office Hadley Centre. (2022). EN4 quality controlled ocean data [Dataset]. https://www.metoffice.gov.uk/hadobs/en4/

Mohrmann, M., Heuzé, C., & Swart, S. (2021). Southern ocean polynyas in CMIP6 models. *The Cryosphere*, 15(9), 4281-4313. https://doi.org/10.5194/tc-15-4281-2021

NEMO Sea Ice Working Group. (2018). Sea ice modelling integrated initiative (SI3) - The NEMO sea ice engine (No. 31). *Zenodo*. https://doi.org/10.5281/zenodo.1471689

Rheinländer, J. W., Smedsrud, L. H., & Nisancioglu, K. H. (2021). Internal ocean dynamics control the long-term evolution of Weddell Sea polynya activity. *Frontiers in Climate*, 3, 718016. https://doi.org/10.3389/fclim.2021.718016

Roquet, F., Madec, G., McDougall, T. J., & Barker, P. M. (2015). Accurate polynomial expressions for the density and specific volume of seawater using the TEOS-10 standard. *Ocean Modelling*, 90, 29-43. https://doi.org/10.1016/j.ocemod.2015.04.002

Schröder, M., & Fahrbach, E. (1999). On the structure and the transport of the eastern Weddell gyre. *Deep Sea Research Part II: Topical Studies in Oceanography*, 46(1-2), 501-527. https://doi.org/10.1016/s0967-0645(98)00112-x

Siebesma, A. P., Bretherton, C. S., Brown, A., Chlond, A., Cuxart, J., Duynkterke, P. G., et al. (2003). A large eddy simulation intercomparison study of shallow cumulus convection. *Journal of the Atmospheric Sciences*, 60(10), 1201-1219. https://doi.org/10.1175/1520-0469(2003)60<1201:alesis>2.0.co;2

Siems, S. T., Huang, Y., & Manton, M. J. (2022). Southern ocean precipitation: Toward a process-level understanding. *Wiley Interdisciplinary Reviews: Climate Change*, 13(6), e800. https://doi.org/10.1002/wcc.800

Sreen, G., Kaleschke, L., & Heygster, G. (2008). Sea ice remote sensing using AMSR-e 89-GHZ channels. *Journal of Geophysical Research*, 113(C2). https://doi.org/10.1029/2005jc003384

Steiger, N., Sallée, J., Ward, B., Azevedo, A., Binase, Z., Ten Doeschat, A., et al. (2022). CTD observation - So-chic cruise 2022 [Dataset]. *SEANOE*. https://doi.org/10.17882/95314

Stoessel, A., Notz, D., Haumann, F. A., Haak, H., Jungclaus, J., & Mikolajewicz, U. (2015). Controlling high-latitude southern ocean convection in climate models. *Ocean Modelling*, 86, 58-75. https://doi.org/10.1016/j.ocemod.2014.11.008

Tsujino, H., Urakawa, S., Nakano, H., Small, R. J., Kim, W. M., Yeager, S. G., et al. (2018). JRA-55 based surface dataset for driving ocean-sea-ice models (JRA55-DO). *Ocean Modelling*, 130, 79-139. https://doi.org/10.1016/j.ocemod.2018.07.002

Wilson, E. A., Riser, S. C., Campbell, E. C., & Wong, A. P. (2019). Winter upper-ocean stability and ice-ocean feedbacks in the sea ice-covered southern ocean. *Journal of Physical Oceanography*, 49(4), 1099-1117. https://doi.org/10.1175/jpo-d-18-0184.1

Zhou, L., Heuzé, C., & Mohrmann, M. (2022). Early winter triggering of the Maud Rise Polynya. *Geophysical Research Letters*, 49(2), e2021GL096246. https://doi.org/10.1029/2021GL096246