Multidecadal Preconditioning of the Maud Rise Polynya Region

. In this paper, we consider Maud Rise polynya formation in a long (250 years) high-resolution (ocean 0.1 ◦ , atmosphere 0.5 ◦ horizontal model resolution) of the Community Earth System Model. We ﬁnd a dominant multidecadal time scale in the occurrence of these Maud Rise polynyas. Analysis of the results leads us to the interpretation that a preferred time scale can be induced by the variability of the Weddell Gyre, previously identiﬁed as the Southern Ocean Mode. The large-scale pattern of heat content variability associated with the Southern Ocean Mode modiﬁes the stratiﬁcation in the Maud Rise re-5 gion and leads to a preferred time scale in convection through preconditioning of the subsurface density, and consequently to polynya formation.

gion and leads to a preferred time scale in convection through preconditioning of the subsurface density, and consequently to polynya formation.

Introduction
A polynya is an open-water area enclosed by sea ice which persists at least for a few months. A famous area for polynya formation is the Maud Rise region in the Weddell Sea. Polynyas occurring at this location are usually referred to Maud Rise Polynyas 10 (MRPs), referring to the bathymetry feature below the ocean water. In 1974, the first MRP was observed and developed into the larger Weddell Polynya in the following years (Carsey, 1980). These polynyas were observed by in situ (Gordon, 1978) and satellite microwave imaging (Carsey, 1980). The mid-1970s polynya was characterised by a sea-ice enclosed open water area (1 -3 × 10 5 km 2 ) located near 0 • E and 65 • S (Gordon, 1978). The interest for polynyas has increased recently because of the appearance of an MRP in the austral winter and spring of 2017, with an open-water area of about 0.5 × 10 5 km 2 (Campbell 15 et al., 2019).
Many observational, modelling and theoretical studies based on the 1974 -1976 polynya have addressed its formation, evolution and decay (Martinson et al., 1981;Parkinson, 1983;Holland, 2001;Gordon et al., 2007;Martin et al., 2013;Latif et al., 2017;Weijer et al., 2017;Kurtakoti et al., 2018). The occurrence of convection in the vicinity of Maud Rise is clearly a generic aspect during MRP formation. Convective events are induced by a destabilisation of the water column and the classical 20 view is that surface salinity anomalies (e.g. through brine rejection) induce convection once the column is preconditioned through subsurface processes (Martinson et al., 1981). In Gordon et al. (2007), it is proposed that a negative phase of the Southern Annular Mode (SAM) causes drier atmospheric conditions (net evaporation) resulting in a more saline surface layer leading to a destabilisation of the water column. Preconditioning can also occur due to the interaction of an ocean eddy and the Maud Rise topographic feature (Holland, 2001) and by stratified Taylor columns (Alverson and Owens, 1996;de Steur et al., The occurrence of Southern Ocean convection is also ubiquitous in many global climate models (GCMs). Martin et al. (2013) show that in the Kiel Climate Model (KCM) a centennial time scale build-up of a heat reservoir at mid-level depths (1000 -3000 m) preconditions the water column affecting convection. Reintges et al. (2017) analysed a variety of (Coupled Model Intercomparison Project phase five, CMIP5) GCMs and found that several models display deep convection events at multidecadal timescales, for example the GFDL-CM (Zanowski et al., 2015;Zhang and Delworth, 2016). Reintges et al. 5 (2017) argue that models with a weak (strong) stable stratification tend to have a shorter (longer) re-occurrence time of deep convection, and demonstrate the importance of sea-ice volume on the average length of both non-convective and convective periods. Dufour et al. (2017) analysed MRP events in two versions (eddy-permitting and eddy-resolving) of the GFDL-CM. The lower resolution version has a weaker vertical background stratification (compared to the high-resolution version) and displays quasi- 10 continuous deep convection events. The effects of eddy-driven heat-and salt transports and the representation of overflows were shown to be crucial for the stratification changes in the Maud Rise region (Dufour et al., 2017). In addition, Weijer et al. (2017) demonstrated in a 100-year model simulation of the Community Earth System Model (CESM) that the MRP is only present in the eddying-version of the same model. The results in Dufour et al. (2017) and Weijer et al. (2017) suggest that ocean eddies (and also overflows) play a significant role in the vertical background stratification (in addition to the effects already mentioned 15 by Reintges et al. (2017)) and therefore affect the time scale of the occurrence of deep convection.
The diverse multidecadal  to centennial (Martin et al., 2013) variability found in climate models gives the impression that convection can occur rather randomly through surface density perturbations once the water column is sufficiently preconditioned. However, recently rather regular multidecadal patterns of internal variability have been found in strongly eddying ocean models in the Southern Ocean, the so-called Southern Ocean Mode (SOM, Le Bars et al. (2016)).

20
Dynamical processes involving eddy mean-flow interaction and eddy-topography interaction (Hogg and Blundell, 2006) are the driving mechanisms of the SOM (Jüling et al., 2018). Recently, van Westen and Dijkstra (2017) have shown that the SOM variability also occurs in a high-resolution ocean-eddying version of the CESM.
In this paper, we pursue the idea of the possible existence of a preferred multidecadal time scale of MRP events in the Southern Ocean through preconditioning and subsequent convection, by analysing the results of an extended CESM simulation 25 as in van Westen and Dijkstra (2017). In this multi-century CESM simulation, several MRP events are found and we analyse the processes involved in these events. In Section 2, information on the CESM simulation and mean properties of the (Antarctic) sea-ice field in the simulation are given. As the Southern Ocean Mode is a relatively recent concept, we provide a brief overview and several characteristics of this mode in Section 3. Analysis of the MRP events in CESM is provided in Section 4. A summary and discussion of the results with the main conclusions are given in the final Section 5. The MRP is defined as the enclosed region within the Antarctic sea-ice pack that has a sea-ice fraction smaller than 15% 10 (Weijer et al., 2017). The observed 2017 MRP and a CESM modelled MRP are both located over Maud Rise (Figure 2b). The MRPs found in CESM turn out to be quite localised in space and to determine that region, we computed a polynya probability density function. First, we define a polynya year to be a year where an MRP appears for at least three months. Secondly, at each MRP, the open-water grid cells (bounded by the 15% sea-ice fraction contour) and sea-ice covered grid cells are given a label of 1 and 0, respectively. This procedure is applied to each polynya year, so non-polynya years are discarded. Taking the 15 average over all the labeled fields (0 or 1), results in the polynya probability density function as shown in Figure 2b. There is a small region near 7.5 • E and 65 • S where there is open water for almost all MRPs. Based on the polynya probability density function, we define the Polynya region (2 • E -11 • E × 66.5 • S -63.5 • S, shown as the dashed outlined region in Figure 2b).
We determined spatial averages over the Polynya region for the temperature, salinity, ocean heat content, sea-ice fraction, sea-ice thickness and the monthly maximum mixed layer depth. The mixed layer depth is defined as the depth at which the interpolated buoyancy gradient matches the maximum buoyancy gradient (this is standard output in CESM, Smith et al. (2010)).

5
The local ocean heat content (OHC), indicated by H k , at model level z k was calculated as: where dz k is the vertical cell length of the model grid at that level. The quantities ρ k , T k and C p,k are the local density, temperature and heat capacity, respectively. The temperature, salinity and pressure dependency are taken into account when calculating the local density and heat capacity (Millero et al., 1980;Sharqawy et al., 2010). The drift in the integrated OHC 10 over the Polynya region (not shown) is smaller compared to that in the time series in Figure 1b.
The propagation of particles into the Polynya region is studied using OceanParcels (Probably A Really Computationally Efficient Lagrangian Simulator, version 2.1.4, Delandmeter and Van Sebille (2019)). In OceanParcels, one can release a set of virtual particles and track their path under the flow using Lagrangian modelling. Particles are passively advected by the time-varying 3D velocity field from the CESM output, either forward or backward in time. When a particle is released, its age 15 is set to zero. We used a time step of ∆t = 1 hour to update the location of each particle (by linear interpolation of the velocity fields, Delandmeter and Van Sebille (2019)) and the output of each particle (i.e. age, depth, latitude and longitude) is stored every 10 days.

The Southern Ocean Mode
In a stand-alone ocean simulation with the POP (not to confuse with the coupled model simulation, i.e. the CESM), Le Bars  Figure 3c (red curve). The first EOF and PC capture 16% of the total variance and the dominant period of the SOM is 40 -50 years (Le Bars et al., 2016;van Westen and Dijkstra, 2017;Jüling et al., 2018).
The variability associated with the SOM can also be measured via the SOM index which is defined as the surface temperature first PC. The lag-correlation at 65 months is r = 0.93 and is significant (95%-confidence level (Cl from now on), taking into account the reduction of the degrees of freedom due to the running mean). The EOF dipole pattern affects the isopycnal slopes in the Southern Ocean (Le Bars et al., 2016) and consequently, via baroclinicity, the Weddell Gyre strength varies with the same period of the SOM (Figure 3e). The first PC and the Weddell Gyre strength time series are 180 • out of phase and have a significant (95%-Cl) correlation of r = −0.90. The SOM is not related to the SAM (Gong and Wang, 1999) as the SAM index is constant in the stand-alone POP (Figure 3e).
Applying the same PCA above for the CESM output, we also find a large-scale pattern of variability in the Southern Ocean ( Figure 3b). The dominant period of variability of the first PC is 25 years (Figure 3d, see also Figure 5b). The SOM index also displays multidecadal variability in the CESM (Figure 3d), but is much more variable compared to that of the stand-alone POP. For example, the lag-correlation between the first PC and the SOM index is much lower in the CESM (r = 0.47, 95%-Cl, 10 36 months) compared to the stand-alone POP (r = 0.93, 95%-Cl, 65 months). The differences in the EOF pattern and SOM index between the stand-alone POP and the CESM are related to atmospheric and sea-ice processes in the CESM. Therefore we propose the SOM * index, which is defined as the upper 1000 m averaged temperature anomaly over the region 30 • W - Figures 3a,b). The SOM * region is based on the maxima of the EOF patterns of both the stand-alone POP and CESM.

15
The multidecadal variability in the CESM is better represented in the SOM * index than in the SOM index ( Figure 3d). For example, from lag-correlation analysis we find that the first PC is highly correlated with the SOM* index (r = 0.80, 95%-Cl, 19 months). For the stand-alone POP, the phase difference between the first PC and the SOM * index is 29 months and r = 0.95 Atmospheric and sea-ice related processes also cause that less variance is captured in the first PC for the CESM (4.1%) compared to the stand-alone POP (16%). The first and second PC are well separated for both the stand-alone POP and CESM (North et al., 1982). Still, the variances are rather low in both models. Low variances in PCA are related to spatial noise in large parts of the ocean, the total number of grid cells and the length of the time series. Part of the spatial noise can be removed 25 when the temperature fields are smoothed (through a 60-month running mean) before applying PCA. In this case, the captured variance by the first PC increases to 28.3% (stand-alone POP) and 15.3% (CESM). Another possibility to increase the variance is by reducing the region (50 in the first EOF (as in Figures 3a, b) when reducing the size of the region or by applying a moving average (not shown). To test whether the variances of PCA of the full field are significantly different from red noise (null hypothesis), we generated red-noise surrogate temperature fields with the same statistical properties as the pre-filtered time series. The variance of the first PC of these surrogate temperature fields is 3.8% and 0.9% at the 99%-Cl (500 surrogate fields) for the stand-alone POP and CESM, respectively. Hence, although the variances of the PCs are relatively small in both the stand-alone POP and CESM, they are significantly different from red noise.
The SAM index (Gong and Wang, 1999)

Results
In subsection 4.1, we first present the characteristics of the MRP events found in the CESM simulation. It turns out that there is a 25 preferred multidecadal variability in MRP formation, possibly linked to the SOM (Le Bars et al., 2016), and this teleconnection is analysed in subsection 4.2. The processes causing the associated convective events, in particular the preconditioning of the density field are analysed in the last subsection 4.3.

Maud Rise Polynyas in the CESM
The September sea-ice fraction (blue curve in Figure 4a Fourier spectra of the mixed layer depth and the first PC (cf. Figure 3d) are shown in Figure 5. The dominant period of convection (i.e. deepening of the mixed layer depth) has a significant (99%-Cl) 25-year period against red noise (Figure 5a). 10 The time series of OHC, temperature, salinity and sea-ice fraction vary with the same multidecadal period as a consequence of convection (spectra not shown). The first PC also varies on a dominant 25-year period which is significant (99%-Cl) against red noise (Figure 5b). The SOM and SOM * indices also display a dominant multidecadal period (20 -30 years, not shown).
The possible connection between the SOM and the mixed layer depth variability is analysed in the following section.

15
As demonstrated in subsection 4.1, both the SOM and the mixed layer depth over the Polynya region display multidecadal variability with the same dominant time scale. In this subsection, the propagation of subsurface heat anomalies towards the Polynya region is analysed. For this analysis, we use the SOM * index instead of the first PC because the first PC also contains variability due to MRP formation. The SOM * region is located far outside the Polynya region.
A lag-correlation analysis between the SOM * index (time series in Figure 3d),H 100 andH 1000 averaged over the Polynya region (time series in Figure 4a) is shown in Figure 6a. There are significant lag-correlations (95%-Cl) between the SOM * index andH 1000 , where the SOM * index leads theH 1000 by 7 years (80 months). At the similar time scale, negative correlations 5 exists betweenH 100 and the SOM * index. The first PC and the SOM index display similar significant lag-correlations with thē H 100 andH 1000 averaged over the Polynya region.
The lag-correlation patterns of the SOM * index and theH 1000 field (Figures 6b and 6c)  year 249. Therefore, we released particles every 30 days in the Polynya region (each spaced by 0.5 • and fixed initial depth, 133 particles per cast) and we backtracked these particles using OceanParcels over the CESM output (model years 150 -250).
The backtracked trajectories demonstrate that the particles propagate along the Weddell Gyre and have an upstream origin.
Rather than showing all the trajectories (as in Figure 7b), we show the particle distribution after 7 years (80 months) of back- The backward propagation of the particles along the Weddell Gyre can also be identified using the colour-coded regions (as in Figures 8a,c,e). Here we determine the fraction of particles inside a specific region at a given age (Figures 8b,d,f). All particles are released in the black-outlined region (this region slightly extends the Polynya region), and hence all particles have 30 zero age. The time-mean flow in the black-outlined region is directed westwards (Figure 7a) which explains the relatively high (≈ 70%) abundance of particles in the blue-outlined region after 2 years of backtracking. Via the blue dashed region, most (> 50%) particles end up in the red region after about 9 -13 years.
The interpretation of these results is that when the SOM is in a positive phase, positive (subsurface) OHC anomalies are found in the SOM * region (Figure 3b). The SOM * index, which is merely a measure of the phase of the SOM, shows significant  correlations with theH 1000 field and also positive correlations are found in the SOM * region (Figure 6b). These subsurface OHC anomalies reach the Polynya region in about 7 years (Figures 6d -f). In this way, the subsurface OHC anomalies are advected along the Weddell Gyre and enter the Polynya region (Figures 7 and 8).

Preconditioning of the Polynya region
The propagation of subsurface OHC anomalies associated with the SOM causes subsurface heat accumulation in the Polynya 5 region and leads to changes in the subsurface heat reservoir. The horizontal advective heat flux, indicated by F k , at model level z k , was calculated as: where u k is the horizontal velocity vector at depth z k . The zonal heat input at the eastern boundary (F East k ) of the Polynya region (i.e. the normal component of F k multiplied by the area perpendicular to the normal) with depth also displays mul- andP Net 1000 ). The cyan curve is the salinity difference between subsurface (S1000) and surface (S100) over the Polynya region. In addition, the blue curve shows the yearly-averaged surface salinity anomaly with respect to model year 210 (S surf ace , magnified by a factor 10) and the dashed blue curve displays the annual maximum sea-ice thickness (Hice, reduced by a factor 2). reservoir.
In model years 210 -215, shortly after the third multiyear MRP event, the heat reservoir is depleted and the temperature difference is relatively small compared to period (model years 225 -230) before the 231 -237 MRP. Between model years 210 -230,F Net 1000 is positive and larger compared toF Net 100 , leading to an increase in the vertical temperature difference. As the 5 temperature difference increases, the total amount of subsurface heat advected out (sum ofF North 1000 ,F South 1000 andF West 1000 ) of the Polynya region also increases, consequently the net heat inputF Net 1000 decreases over time. The time series (model years 150 -250) of the temperature difference andF Net 1000 are significantly anti-correlated (95%-Cl, r = −0.88). The build-up of the subsurface heat reservoir weakens the stratification over the Polynya region. During polynya formation, heat is vertically mixed causing a near zero temperature difference which indicates that the upper 1000 m is well mixed. Once convection ceases, the  (Campbell et al., 2019). However, there are observations near Maud Rise indicating subsurface heat accumulation over longer periods of time (Smedsrud, 2005;Fahrbach et al., 2004Fahrbach et al., , 2011.
Preconditioning of the Polynya region also occurs by advection of salinity anomalies into the Polynya region and the advec-15 tive salt fluxes are given by: The salt input at the eastern boundary (P East The eastern boundary contributes to subsurface salt accumulation in the Polynya region (Figure 9b), the other three boundaries (not shown) advect salt out of the Polynya region. The sign ofP Net 100 andP Net 1000 are opposite, leading to an increase in the salinity difference between the subsurface (S 1000 , 200 -1000 m) and surface (S 100 , upper 100 m) over the Polynya region prior to MRP formation (Figure 9d). An increase in the vertical salinity gradient, strengthens the stratification of the Polynya region.

25
A saline surface layer (due to brine rejection) weakens the stratification, potentially leading to convection (Martinson et al., 1981). Therefore, we determined the surface salinity and sea-ice thickness averaged over the Polynya region which are also shown in Figure 9d. The yearly-averaged surface salinity time series is expressed as an anomaly with respect to model year 210.
For the sea-ice thickness time series, we retained the annual maximum. Prior to MRP formation, the annual maximum sea-ice thickness is increasing over time leading to more brine rejection but this does not result in any increase of the surface salinity 30 over the same period. Taking the average over the months April -July (i.e. sea-ice growth season) for the surface salinity also results in a decrease of surface salinity values prior to MRP formation. Changes in the vertical salinity gradient are mainly determined by horizontal (sub)surface advection of salt. Stratified Taylor columns over Maud Rise also contribute to preconditioning of the region (Alverson and Owens, 1996;de Steur et al., 2007). The water column over Maud Rise is weakly stratified compared to the surroundings but there is no (deep) convection prior to MRP formation (Figure 10a). Here we determined the potential density (standard output of the CESM) and mixed layer depth between 65 • S -64 • S and across 20 • W -30 • E; this section crosses the Maud Rise summit (near 3 • E). There are also Taylor columns over Maud Rise during other years, as shown by the potential density at 900 m depth in Figure 10b.

5
During the polynya years, the potential density at 900 m depth (as well as for other depth levels) tend to zero over Maud Rise indicating that the water column overturns. Gordon et al. (2007) suggest that a negative phase of the SAM contributes to the preconditioning of the region, however, the SAM is (strongly) positive before the first (158 -159), third (205 -209) and fourth (231 -237) multiyear MRP event in the CESM (Figure 3f).
Based on this analysis of the CESM results, the preconditioning of the density field seems to be most affected by the 10 subsurface heat anomalies and Taylor columns over Maud Rise. Positive subsurface heat and salt anomalies, advected via the Weddell Gyre, increase the vertical gradient in temperature and salinity over the Polynya region, respectively. The increased vertical temperature (salinity) gradient weakens (strengthens) the stratification near Maud Rise.

Summary and Discussion
We analysed the last 100 years of model output from a 250-year control simulation of a high-resolution version of CESM under a repeated seasonal forcing of the year 2000. We find four multiyear MRP events and, in each event, (deep) convection causes vertical mixing of anomalous subsurface heat towards the surface where it melts the sea ice and leads to the formation of the MRP. These processes of the formation and lifecycle of the MRP are broadly in agreement with the classical view as in Martinson et al. (1981) and those described from other model results in Martin et al. (2013), Dufour et al. (2017) andReintges et al. (2017). In most model studies deep convection occurs rather randomly while in the CESM deep convection occurs about every 25 years.
In the CESM, we find that the multidecadal preconditioning of the subsurface density can be linked to an intrinsic dynamical 5 ocean mode in the Southern Ocean, the SOM, which is dominantly caused by eddy-mean flow interaction (Jüling et al., 2018).
A positive phase of the SOM leads to positive subsurface OHC anomalies in the South Atlantic Ocean (Le Bars et al., 2016).
These anomalies enter the Weddell Gyre near 30 • E and eventually reach the Polynya region, where the anomalies cause deep convection and sea-ice melt. Hence, the frequency of occurrence of the MRP events is related to the SOM variability through the advection of the subsurface OHC anomalies in the Weddell Sea. In this case, a preferred frequency of convective events is 10 induced through preconditioning, which is around 25 years in the CESM.
Although preconditioning of the density field in the Maud Rise region is connected to subsurface advection of heat and salt anomalies and Taylor columns, it does not strictly imply polynya formation. However, it will favour the occurrence of MRPs due to atmospheric variability, such as intense winter storms. For example, the formation of the 2016 -2017 MRP was initiated by an increased storm frequency under the influence of a positive SAM index (Campbell et al., 2019). In the CESM, the SAM 15 index is in a positive phase prior to MRP formation, but the initiation process of the MRPs is out of the scope of this study.
Observations are unfortunately too sparse to falsify the hypothesis of a preferred multidecadal variability in the Southern Ocean. Reanalysis of sea surface temperatures are available to derive a historical SOM index, but the time series is too short to adequately detect the SOM. Relatively low sea-ice fractions were observed over Maud Rise in the mid-1970s, 1980, 1994and 2016(Comiso and Gordon, 1987Lindsay et al., 2004;Gordon et al., 2007;Campbell et al., 2019), suggesting a 20 15 -20 year return period of MRP formation which is a bit shorter compared to the 25-year period of the CESM MRPs. The background state of the Southern Ocean sets the re-occurrence time for deep convection  which explains the differences in the return period of MRP formation. Besides, the stand-alone POP has a different background state compared to that of the CESM, leading to much longer re-occurrence time of deep convection in the stand-alone POP compared to the CESM (Jüling et al., 2018). Nonetheless, if the view is correct that a preferred multidecadal frequency (∼20 years) in the 25 Southern Ocean exists, it would imply that a next MRP event, if not affected by climate change (De Lavergne et al., 2014), can be expected before the mid-2030s.
Acknowledgements. The authors thank Michael Kliphuis (IMAU, UU) for performing the CESM simulations. The computations were performed on the Cartesius at SURFsara in Amsterdam, the Netherlands. Use of the Cartesius computing facilities was sponsored by the Netherlands Organization for Scientific Research (NWO) under the project 15552. The data from the model simulation used in this work 30 is available upon request from the authors. The NOAA/NSIDC provided the satellite sea-ice products (http://nsidc.org/data/G02202). The OceanParcels framework can be obtained online (oceanparcels.org) and we thank Peter Nooteboom (IMAU, UU) for his assistance in setting up this code for the CESM data.