Multidecadal Polynya Formation in a Conceptual (Box) Model

Maud Rise Polynyas (MRPs) form due to deep convection, which is caused by static instabilities of the water column. Recent studies with the Community Earth System Model (CESM) have indicated that a multidecadal varying heat accumulation in the subsurface layer occurs prior to MRP formation due to the heat transport over the Weddell gyre. In this study, a conceptual MRP box model, forced with CESM data, is used to investigate the role of this subsurface heat accumulation in MRP formation. Cases excluding and including multidecadal varying subsurface heat and salt fluxes are considered and multiple polynya events 5 are only simulated in the cases where subsurface fluxes are included. The dominant frequency for MRP events in these results, approximately the frequency of the subsurface heat and salt accumulation, is still visible in cases where white noise is added to the freshwater flux. This indicates the importance and dominance of the subsurface heat accumulation in MRP formation.

the background stratification and hence precondition the water column. These processes also cause a general halo of relatively low sea-ice concentrations over Maud Rise (Lindsay et al., 2004).
Climate models provide the opportunity to study deep convection and consequently MRP formation. From several climate models, it is known that deep convection in the Southern Ocean varies on multidecadal to multicentennial time scales (Martin et al., 2013;Zanowski et al., 2015;Latif et al., 2017;Weijer et al., 2017). Several climate models show subsurface heat 5 accumulation prior to deep convection, e.g. in the Kiel Climate Model (KCM) (Martin et al., 2013), the Community Earth System Model (CESM) (van Westen and Dijkstra, 2020a) and the Geophysical Fluid Dynamics Laboratory Climate Model (GFDL CM2-0) (Dufour et al., 2017). Through buoyancy gain in the subsurface layer, deep convection is induced, which results in MRP formation (Martin et al., 2013;Latif et al., 2017;Reintges et al., 2017). In the KCM, for example, a stronger stratification results in a longer period for deep convection, because more buoyancy gain is necessary to overcome the more 10 stable stratification Reintges et al., 2017). Another important feature is model resolution as shown by Weijer et al. (2017): MRPs were found in a high-resolution (0.1 • ) version of the CESM, whereas in the low-resolution (1 • ) version of the same model no MRPs were simulated. Dufour et al. (2017) used the GFDL CM2-0 model with a nominal ocean grid spacing of 0.25 • and 0.1 • and they show that the occurrence of deep convection itself is not sufficient to create MRPs. If the subsurface heat reservoir cannot supply enough heat to melt all the sea ice, an MRP will not form.

10
The MRP box model consists of two vertically stacked boxes with a constant depth and within each box the ocean properties (e.g. temperature and salinity) are uniform. The model simulates the development of temperature (T), salinity (S), and sea-ice thickness (δ) in each box under surface and subsurface forcing. The depth of the entire water column is H with the surface layer having a depth of h 1 and the subsurface layer a depth of h 2 .
The model has four different flow regimes, below referred to as regimes, which are differentiated on sea-ice cover (sea-ice 15 free versus sea-ice covered) and static stability (two layered versus mixed). Whenever stable/unstable is mentioned below, we refer to static stability of the water column, so no dynamical instabilities. There are the sea-ice free regimes I and II, and the sea-ice covered regimes III and IV. Regimes II and IV are stably stratified (ρ 1 < ρ 2 ), and regimes I and III are mixed with one uniform density over the entire depth (Fig. 1). The subscripts 1 and 2 correspond to the surface and subsurface layer, respectively.

20
Over time the model state may transit through these four regimes under the influence of (seasonal) forcing. The four different regime transitions are indicated by the arrows in Fig. 1, i.e., (a) From sea-ice covered regimes to sea-ice free regimes due to complete melt of the sea ice (δ = 0) (regime IV → II or regime III → I), (b) from sea-ice free regimes to sea-ice covered regimes, because the surface layer reaches freezing temperature and sea ice 25 starts to form (T 1 = T f ) (regime II → IV, and regime I → III), (c) from stable, two layered regimes to unstable, mixed regimes, because the density of the surface layer is equal to or larger than that of the subsurface layer (ρ 1 ≥ ρ 2 ) (regime II → I, and regime IV → III). The water column becomes unstable and mixes through overturning. Temperature and salinity are uniform over the entire layer and indicated by T and S instead of T n and S n , and 30 (d) from unstable, mixed regimes to stable, two layered regimes, because of stabilisation of the water column due to a decreasing density of the mixed layer (regime I → II, and regime III → IV). The precise conditions for the regime transitions are shown in Appendix A. It should be noted that for the model to switch between regimes I (mixed, sea-ice free) and III (mixed, sea-ice covered) the entire water column should reach freezing temperature, which is physically not realistic. Therefore, this transition does not exist in the model (regime I regime III) The model is forced at the surface by a freshwater flux F , and by a heat flux Q ia for sea-ice covered regimes and Q oa for open-ocean regimes. ::: The ::::::::: freshwater :::: flux :: F : is :::::::: modelled ::: as : a :::::: virtual ::: salt :::: flux :: by :::::::::: multiplying : it :::: with :: a :::: base :::::: salinity ::: S 0 . : Both the F T 2 and F S2 for the subsurface layer) which depend on a background value (T b1 and S b1 for the surface layer, and T b2 and S b2 for the subsurface layer) and a relaxation timescale (τ ).
The heat and salt transfer between the layers are modelled using exchange coefficients (K T and K S ), which account for upwelling, turbulent exchange and diffusion. In sea-ice covered regimes there is a heat flux present between the sea ice and the underlying layer. This heat flux is modelled using a turbulent exchange coefficient (K). Brine is rejected during sea-ice growth 5 and freshwater is added to the surface layer during sea-ice melt, brine rejection and sea-ice melt are modelled using a constant representing the salinity difference between sea ice and seawater (σ) and the rate of sea-ice growth (dδ/dt). The density for each layer is determined from a simple linear equation of state: where the constants α and β are the thermal expansion and haline contraction coefficients, respectively.

10
The governing equations for each regime are then: Regime I: Regime II: Regime III: :::: Regime IV: :::::::::::: In these equations, C p is the specific heat of seawater with a reference density ρ 0 . Sea ice has a density of ρ i and a latent 15 heat of melting/freezing indicated by L. The standard values of all parameters used in the model are presented in Section 2.3.
The ODE15s solver is a variable-step, variable-order solver based on an algorithm by Klopfenstein (1971) using numerical differentiation formulas (NDFs) orders 1 to 5. Tolerances for the absolute and relative error are used to increase the accuracy of the model; these tolerances are set to 10 −10 and 10 −8 , respectively.

MRP Box Model Setup and Case Description
The original polynya model of Martinson et al. (1981) is obtained from the model formulation above by setting the horizontal advective fluxes (F T 1 , F S1 , F T 2 and F S2 ) to zero and by setting the subsurface layer values of temperature and salinity to 25 constant values. We consider two cases of this configuration of the model, which only differ by the value of K S used. The higher value of K S (case MKH: Martinson K S high) results in more salt transfer from the subsurface layer to the surface layer (than in case MKL (Martinson K S low) with the lower value of K S ) increasing the density of the surface layer, making it more susceptible to overturning (Table 1).
The three cases with a dynamic subsurface layer are differentiated on the inclusion of the different components of the Parameter values for each case are displayed in Table 2. The parameter values are either taken from Martinson et al. (1981), or based on the CESM simulation of van Westen and Dijkstra (2020a), or they are determined through tuning of the model.

5
For the CESM simulation, we determined spatial-averaged quantities over the Polynya region (2 , where an MRP forms in the CESM (van Westen and Dijkstra, 2020a). The aim of the model is to investigate multiple polynya events, and therefore it is necessary to tune the stratification. The stratification can be either too strong and no overturning occurs, or the stratification is too weak, and the water column overturns each year. To obtain multiple polynya events the heat and salt fluxes between the two layers and between the sea ice and the surface layer are tuned.

10
The typical depth of the layers has been determined from the CESM simulation of van Westen and Dijkstra (2020a). The depth of the surface layer (h 1 ) is set to 160 m, because the potential density shows a clear homogeneous layer below 160 m   Figure 2).  The turbulent exchange coefficient K, and the exchange coefficients K T and K S have been used as tuning parameters for the different cases. The coefficient K is set to 1 × 10 −4 ms −1 for all cases (in Martinson et al. (1981) this value is 3 × 10 −4 ms −1 ).

10
The values of K T and K S per case are shown in Table 1. To compare the magnitude of these parameters with values used in literature the values need to be converted from ms −1 to m 2 s −1 , which is the usual unit for vertical diffusivity parameters. This is done by multiplying these values with the depth of the surface layer (i.e. 160 m). This results in values between 2.2 × 10 −4 m 2 s −1 and 5 × 10 −4 m 2 s −1 . Comparable values are found in a model study of Dufour et al. (2017) for this same location and in observations (Shaw and Stanton, 2014). The values used in this study are up to a factor 10 larger than the values 15 used in Martinson et al. (1981) (K T = 7 × 10 −7 ms −1 and K S = 10 −7 ms −1 ).
Initial conditions of the model affect the long-term behaviour of the model. Besides, the initial regime of the model should match with initial conditions. For example, when starting in a sea-ice covered regime, the initial conditions for sea ice should be δ > 0. Another important initial condition is the initial stratification. If the initial stratification is too weak, the model will overturn each year. In this specific case, it is not possible to study multidecadal variability in polynya formation. Therefore, 20 each model simulation is initiated on January 1 with the following conditions: T 1 = 0.1 • C, S 1 = 34.2 g/kg and δ = 1 m.

Forcing Conditions
The MRP box model is forced by a monthly varying heat flux (Q oa or Q ia ) and a monthly varying freshwater flux (F ) that are repeated for each model year (Fig. 3c,d, Table 3). The monthly-averaged heat fluxes are retained from the CESM. The fluxes are spatially averaged over the polynya region (2 • E -11 • E × 63.5 • S -66.5 • S, black outlined region in Fig. 2) defined in van Westen and Dijkstra (2020a). The surface heat flux strongly increases when the sea-ice fraction is lower than 0.5 in 10 the CESM (not shown). Therefore, sea-ice fractions lower (higher) than 0.5 represent a sea-ice free (covered) regime with the ocean-atmosphere (sea ice-atmosphere) heat flux noted by Q oa (Q ia ). Note that the MRP box model has a discrete sea-ice fraction of either 0 or 1. The monthly averaged heat fluxes are linearly interpolated in time. We do not use the sea ice-ocean flux from CESM, since this flux is already represented in the model equations (e.g. by the term K(T − T f ) in equation (4)).
There is a tendency of more evaporation (decrease in F ) when an MRP forms. Therefore, the model uses different values  Values for the four horizontal advective fluxes (F T 1 , F S1 , F T 2 and F S2 ), background temperatures (T b1 and T b2 ) and salinities (S b1 and S b2 ) are obtained from the CESM simulation. The first layer uses a constant background temperature (T b1 = −0.33 • C). The constant background salinity (S b1 ) for the surface layer is slightly changed to tune the model (from 34.5 to 34.4818 g/kg). This is necessary to be able to simulate multiple polynya events. In the CESM, the background temperature and salinity of the subsurface layer (200 -1000 m depths) are periodically varying as shown in Fig. 4. The dominant period in 5 CESM is 25 years (not shown). This translates into a prescribed 25-year cycle for F T 2 and F S2 . The time mean of T b2 and S b2 are also shown in Fig. 4. ::: The ::::: effect :: of :::::::: stratified ::::: Taylor :::::::: columns :: are :::::::::: assimilated :: in ::: the ::::::::: subsurface :::::::::: temperature ::: and ::::::: salinity :::: time ::::: series. : These horizontal fluxes are used to make sure the water masses in the box do not drift away from the surrounding water masses. We have used fitted background states since we are using a highly idealised model incapable of reproducing the

Yearly Repeated Cycles
The MRP box model displays three types of yearly cycles, which are shown schematically in Fig. 5 for the surface box temperature and salinity. Note that the salinity axis has no values. If actual values would have been used, the cycles would overlap (see e.g. in Fig. 7c). In Fig. 5, each cycle starts at 'A' and follows in alphabetical order where each letter stands for a regime transition. We use the following definition for MRP formation in this box model: An MRP has formed when the sea ice 5 has melted away completely while there is still ocean cooling (negative heat flux).
For the 0-overturn cycle, the model cycles between regime II (stable, sea-ice free) and IV (stable, sea-ice covered). At 'A', the model transits from regime II → IV because the surface layer reaches the freezing temperature. During sea-ice growth, brine is rejected resulting in an increase of the surface salinity. During austral spring, the sea ice melts leading to an increase of the freshwater flux and consequently salinity levels decrease. At 'B', the model transits back to regime II since all the sea 10 ice has melted. Next, the surface layer temperature T 1 is mainly controlled by the atmospheric heat flux, which is positive (negative) during austral summer (autumn). When the model is cooled down to freezing temperatures, the model transits in 'A' to regime IV and the (seasonal) cycle continues.
For the 1-overturn cycle, where the model overturns once, the surface layer also reaches freezing temperatures in 'A' and brine is rejected. The difference here is that the water column becomes statically unstable (ρ 1 ≥ ρ 2 ) in 'B' and enters the 15 overturning regime III (mixed, sea-ice covered) for a relatively short period (in the order of minutes) after which in transits back to regime IV in 'C'. Due to vertical mixing, relatively warm and saline subsurface water is mixed upwards and melts the sea ice; a polynya forms. During polynya formation, the model transits to regime II in 'D'. Similar to the 0-overturn cycle, the surface layer is controlled by the atmospheric heat flux for the remaining part of the year (austral winter -autumn). This cycle is similar to the one reported in Martinson et al. (1981).
For the 2-overturn cycle, the model transits through all four regimes. The first part of the cycle (A -D) is similar to the 1overturn cycle and a polynya forms in 'D' but considerably less sea ice forms between 'A' and 'B' compared to the 1-overturn cycle. After 'D', the surface layer is cooled during austral winter, but the difference here is that the surface layer becomes static 5 unstable in 'E' and switches to regime I (mixed, sea-ice free). After the model state becomes statically stable in 'F', it stays stable for the remaining part of the year.

Cases MKL and MKH
When the MRP box model was configured with parameter values and numerical schemes as reported in Martinson et al. (1981), their results could not be reproduced. Unfortunately, there is an incomplete parameter documentation in Martinson et al. (1981).
parameters are slightly different than the ones in Martinson et al. (1981)), the MRP box model is spun-up :::: spun ::: up for 75 years and continued for another 25 years (model years 76 -100).
The results for the cases MKL and MKH are shown in Fig. 6a-c and Fig. 6d-f, respectively. For both cases, the density of both layers, the sea-ice thickness and a T-S diagram are shown for the last 25 years. MKL remains in the 0-overturn cycle ( Fig. 5) while in MKH (larger K S ), a 2-overturn cycle is found. In Fig. 6b and e, the sea-ice thickness shows the same (yearly) 5 cycle every year. For temperature and salinity in the surface layer this can be seen in Fig. 6c and f indicating that there are no variations with a period larger than 1 year. In case MKH, there is a little sea-ice growth each year followed by overturning and subsequent sea-ice melt. Using our definition of polynya formation, we find that a polynya forms each year, and therefore case MKH can be considered as having one long polynya period. In summary, for both cases MKL and MKH only yearly cycle solutions are found under the given forcing. Both solutions do not correspond with MRP events in observations, since no 10 persistent MRP is found over a few years (Campbell et al., 2019).
In the CESM results of van Westen and Dijkstra (2020a) Fig. 8b). In a 25-year cycle, the water column overturns when the subsurface density is approaching its minimum. We can also see that years with overturning can be separated by a year without overturning (e.g. in Fig. 8 in year 76 and 78 there is overturning, but not in year 77). The 30 subsurface processes also influence the characteristics of the surface layer (Fig. 8c). In cases MKL and MKH the yearly cycles overlap each other but in the case PFB the yearly cycles are different as a response to the subsurface heat and salt accumulation which is also seen in the CESM simulation (Fig. 7). Whereas PFB uses both subsurface fluxes (heat and salt), case PFH ( Fig. 9) only uses a time-varying subsurface heat flux forcing. By thermal expansion, the subsurface density decreases and the water column becomes statically unstable in model year 86. During vertical mixing, relatively low sea-ice fractions are found compared to the stable model years. Of the simulated 25 years, there are 13 non-polynya years and 12 polynya years. Note that here we do have one long polynya period, which is different from what we see in the PFB case (Fig. 8). In case PFS (Fig. 10), only a time-varying salt subsurface flux forcing is considered. By haline contraction, the subsurface density increases and the water column is statically stable during relatively high levels of subsurface salinity. When subsurface salinity levels decrease over time, the water column becomes unstable and In all three cases with subsurface forcing, the MRP box model is able to simulate the general features also seen in the CESM simulation (van Westen and Dijkstra, 2020a). These cases show a repeating 25-year cycle, which is the same period as the period of the subsurface forcing and the same period as seen in CESM. Where CESM has more non-polynya years 10 than polynya years case PFS have more polynya years, and case PFH has as many non-polynya years as polynya years. Case PFB, however, approaches the ratio non-polynya years versus polynya years seen in the CESM. Besides this difference, also the timing of the first overturn after a non-polynya period is different with respect to CESM. In CESM the first overturn occurs approximately 6 years after the subsurface heat and salt accumulation have reached their maximum. PFS overturns in approximately the same year, while PFB overturns 3 years later. Case PFH differs most, since it overturns 8 years earlier, and even before the subsurface heat accumulation has reached its maximum. These differences are probably caused by the idealizations in the MRP box model, and most likely due to the representation of mixing in this model, compared to that in CESM.

5
In this section we analyse whether the multidecadal MRP variability, as found for the cases PFB, PFH and PFS, is robust under the influence of atmospheric variability, such as intense winter storms. This atmospheric variability is incorporated into the MRP box model by adding white noise to the surface freshwater flux. White noise was added as in Eq. (6): where F N is the freshwater flux with noise, F is the freshwater input without noise as in Fig. 3d, σ N the standard deviation of the noise (0.6613 mm/day), and r is a random draw from a standard normal distribution on every time step. The standard deviation of the noise was determined from the CESM simulation. Figure 11 displays the spectral power of the variables T 1 (Fig. 11a), S 1 (Fig. 11b), and δ (Fig. 11c) for case PFB. For all variables, the percentile, mean and median, the dominant period is about 25 years, the same period as the subsurface forcing 5 and is clearly visible in all variables (Fig. 11a). The single ensemble member also shows a dominant period of approximately 10 years, showing that the noise can also induce shorter periods of convection. Using the same white noise forcing in case MKL yields no dominant multidecadal period (not shown). As was seen in Section 3.3, the MKL case remains in 0-overturn cycle without polynyas (Fig. 6). Atmospheric noise can cause polynyas in the MKL case. However, when MKL is forced into a polynya state, it cannot be forced out of the polynya state; a polynya forms every year. This means MKL with noise will in time show the same behaviour as MKH without noise. The combination of PFB and MKL shows that atmospheric variability can alter the timing of polynya formation but the dominant period is set by subsurface fluxes of heat and salt.

Summary and discussion
In this study, an extended version of the idealized box model of Martinson et al. (1981), was used to investigate the importance of surface forcing and subsurface forcing on Maud Rise Polynya (MRP) formation, with the aim to understand the CESM original paper), the qualitative behaviour of the MRP box model (reduced to the case used in Martinson et al. (1981)) was the same.
The results for the cases MKL and MKH (close to the case in Martinson et al. (1981)) show that deep convection is caused by brine rejection in a preconditioned surface layer. Brine rejection causes a rapid increase of density in the surface layer. The results ( Fig. 6d-f) clearly show that this eventually induces deep convection. However, brine rejection alone cannot explain 5 observed multiple polynya events (e.g. the 1970s and 2017 events), since brine rejection is present in all years with sea-ice growth, and not all years show deep convection and subsequent polynya formation (Fig. 6a-c). In MKL, there is no overturn (Fig. 6a-c) and in case MKH, two overturns occur each year (Fig. 6d-f). In Martinson et al. (1981), it was also shown that in time the model state will reach a yearly repeating cycle, either in regimes II (sea-ice free, stable) and IV (sea-ice covered, stable) (0-overturn cycle), or in regimes I (sea-ice free, mixed) and II (sea-ice free, stable) (2-overturn cycle). The 0-overturning the surface layer is too weak to overcome the stratification. In the second solution (two overturns each year) this salt transfer towards the surface is increased, which causes static instability in the water column with mixing as a result.
In climate models, such as CESM, the water column stabilises after deep convection when the heat and salt reservoirs are depleted. This depletion leads to stabilisation of the water column by increasing the density of the subsurface layer through heat depletion. This physical process is missing in cases MKL and MKH, and therefore the MRP box model state is not able 5 return to a non-polynya regime. However, in Martinson et al. (1981), a solution was shown with overturns in the first years, after which a stable, non-polynya cycle appeared. These overturns in the first years are the result of the initial conditions. When the model is forced by the CESM surface forcing, either a polynya forms each year (MKH) or the model stays in a stable nonpolynya state (MKL). Clearly, the cases MKL and MKH cannot capture the behaviour of the CESM as found in van Westen and Dijkstra (2020a).

10
In cases PFB (Fig. 8), PFH (Fig. 9) and PFS (Fig. 10), subsurface forcing derived from CESM was prescribed in the MRP box model. The results showed periodic polynya events with the same dominant period as seen in van Westen and Dijkstra (2020a) caused by the periodic subsurface forcing. The subsurface forcing preconditions both the subsurface and the surface layer, after which brine rejection is essential to induce deep convection. Note that in the CESM results of van Westen and Dijkstra (2020b), brine rejection was less important as convection was initiated at the subsurface. Because this subsurface-15 iniated convection in CESM is spatially localized, we probably do not find it in the MRP box model. The finding of van Westen and Dijkstra (2020b) that the subsurface heat flux is more dominant than the subsurface salt flux is not confirmed in the box model, since all cases PFB, PFH, and PFS show comparable behaviour. The subsurface heat flux, however, is expected to be dominant as it influences every quantity (T 2 , ρ 2 , T 1 , ρ 1 , δ, S 1 , S 2 ) in the model. The subsurface salt flux only affects the density and salinity of both layers, and therefore is expected to have a much smaller influence on the results. Kaufman et al. 20 (2020) also found heat build up in the ocean and attributed this to reduced heat loss under sea-ice covered conditions. Ocean heat advection actually seemed to counteract the heat build up. In our model, such a situation does not occur, since T 2 is always larger than T b2 because the subsurface layer is losing heat to the surface layer. The MRP box model is able to capture the general features of MRP formation as seen in van Westen and Dijkstra (2020a) and shows the importance of the subsurface forcing. However, the model is still too idealised to accurately capture the precise 25 MRP formation processes in the CESM simulation. The asymmetry in the non-polynya regime versus the polynya regime was poorly captured is cases PFH and PFS. The asymmetry in case PFB compared best to the CESM but still showed relatively more polynya years compared to the CESM simulation. This is probably due to the difference in how vertical mixing is represented.
In the MRP box model the layers are either in a stably stratified configuration with a constant layer depth, or they are completely mixed. In van Westen and Dijkstra (2020a), a KPP boundary mixed layer scheme is used. Representing the growth of the mixed 30 layer more accurately would improve the model, and possibly would lead to a better representation of the asymmetry between the two regimes. When the mixed layer is allowed to grow more gradually, a lag is introduced in the system. This will delay the formation of a polynya. Due to this instant mixing, both temperature and salinity in the surface layer increase instantly.
This results in large differences after overturning between the MRP box model results and the CESM simulation results (van Westen and Dijkstra, 2020a).

35
they suggested that surface processes are responsible for polynya formation, a view that is still widely supported nowadays.
What we have shown is that their model ::::::: adjusted :::::: model :: to :: the :::::: Maud :::: Rise ::::: region : is not capable of simulating multiple polynya events as seen in observations. When the model is extended with variable subsurface heat and salt accumulation, the model is capable of simulating multiple MRP events and able to qualitatively reproduce the CESM results of van Westen and Dijkstra 5 (2020a). Our study suggests that surface related processes cannot completely explain MRP formation nor long-term MRP variability, and that subsurface advective processes need to be taken into account.