## Abstract

The average predictability time (APT) method is used to identify the most predictable components of decadal sea surface temperature (SST) variations over the Southern Ocean (SO) in a 4000-yr unforced control run of the GFDL CM2.1 model. The most predictable component shows significant predictive skill for periods as long as 20 years. The physical pattern of this variability has a uniform sign of SST anomalies over the SO, with maximum values over the Amundsen–Bellingshausen–Weddell Seas. Spectral analysis of the associated APT time series shows a broad peak on time scales of 70–120 years. This most predictable pattern is closely related to the mature phase of a mode of internal variability in the SO that is associated with fluctuations of deep ocean convection. The second most predictable component of SO SST is characterized by a dipole structure, with SST anomalies of one sign over the Weddell Sea and SST anomalies of the opposite sign over the Amundsen–Bellingshausen Seas. This component has significant predictive skill for periods as long as 6 years. This dipole mode is associated with a transition between phases of the dominant pattern of SO internal variability. The long time scales associated with variations in SO deep convection provide the source of the predictive skill of SO SST on decadal scales. These analyses suggest that if the SO deep convection in a numerical forecast model could be adequately initialized, the future evolution of SO SST and its associated climate impacts are potentially predictable.

## 1. Introduction

Over the past decade, the observed sea surface temperature (SST) in the Southern Ocean (SO) did not increase (e.g., Latif et al. 2013; Zhang et al. 2017b) but instead exhibited cooling anomalies. The associated Antarctic sea ice extent showed an expansion, with a record maximum occurring in September 2012 (e.g., Cavalieri and Parkinson 2008; Comiso and Nishio 2008). In the meantime, the Southern Ocean subsurface (below 500 m) warmed considerably (Purkey and Johnson 2010, 2012). The slowdown in the rate of SO warming cannot be attributed to a decrease in greenhouse gas emission from human activity. Climate models forced by observed temporally varying radiative forcing do not reproduce the observed cooling around the Antarctic, but instead simulate a slow but steady warming and Antarctic sea ice loss (Purich et al. 2016). It is therefore likely that internal variability is contributing to the declining SSTs in the SO (Cane 2010; Zunz et al. 2013; Polvani and Smith 2013). However, the extent to which such SO internal climate variability can be simulated and hence predicted on decadal time scales is still not known.

Decadal predictions are in high demand by decision makers who help plan infrastructure investments and resource rearrangements (Cane 2010). The scientific basis of decadal prediction should be built firmly before this demand can be met. A first step is to estimate whether there is a potentially predictable component on decadal scales. Decadal predictability is commonly estimated by two approaches: prognostic and diagnostic approaches (e.g., Pohlmann et al. 2004). In the prognostic approach, decadal predictability is evaluated based on an atmosphere–ocean fully coupled model (AOGCM) initialized by identical oceanic and perturbed atmospheric conditions. The spread within the ensemble is interpreted as an estimate of predictability. Previous studies further extended this method to decadal hindcasts/forecasts that are initialized with observations. Pioneering studies using prognostic method primarily focused on the North Atlantic and North Pacific Oceans where the observations are more numerous and can better be used for model initialization (Keenlyside et al. 2008; Smith et al. 2007; Robson et al. 2012; Yeager et al. 2012; Yang et al. 2013; Msadek et al. 2014; Mochizuki et al. 2010; Meehl and Teng 2012). These model results suggested that the observation-based initial conditions improve skill in the North Atlantic and, to a lesser extent, the North Pacific.

Compared to prognostic approaches, the diagnostic approaches are easier to carry out since they do not require extensive data for initializing prediction models. Diagnostic predictability can be evaluated by various statistical methods, including examining eigenmodes of a linear inverse model (LIM) (Newman 2007), examining the growth of optimal perturbations (Zanna et al. 2012), and investigating the potential predictability variance fraction (ppvf) (Boer 2004; 2011). These statistical tools can identify where and on what time scale the variables have potential high predictability. Using multimodel ensemble data participating in the Coupled Model Intercomparison Project (CMIP), Boer (2004) found that the largest potential predictability on decadal scales is predominately over the high-latitude oceans, particularly in the SO and North Atlantic. These diagnostic approaches might serve as a useful benchmark for decadal predictions that are based on observation-initialized numerical models.

Given the dearth of long-term observations over the SO, we choose to use a diagnostic method to investigate the potential predictability of decadal-scale SO SST variations by taking advantage of a long control integration of the GFDL CM2.1 model. The method we used here is called average predictability time (APT), as proposed by DelSole and Tippett (2009a,b). The APT method finds the most predictable patterns. One advantage of this technique is that it can capture predictable features that contribute little to total variance growth or cannot be expressed as oscillatory modes (DelSole et al. 2013). The main goal of the current study is to examine the leading predictable components of SO SST and the associated climate impacts within a long control simulation of the GFDL CM2.1 model. The physical processes contributing to this predictability are also investigated. We hope our diagnostic analysis can provide a useful reference point for future SO decadal forecasts using observation-initialized numerical models.

## 2. Models and methods

### a. Coupled model

The long-time integrated control run we used in the present paper comes from the Geophysical Fluid Dynamics Laboratory (GFDL) Coupled Model version 2.1 (CM2.1; Delworth et al. 2006). The CM2.1 model has an atmospheric horizontal resolution of 2° × 2°, with 24 levels in the vertical. The ocean and ice models have a horizontal resolution of 1° in the extratropics, with meridional grid linearly decreasing to ⅓° near the equator. The ocean model has 50 levels in the vertical, with 22 evenly spaced levels over the top 220 m. A 4000-yr control simulation is conducted with atmospheric constituents and external forcing held constant at 1860 conditions. We perform analyses using the last 3000 years of the simulation (1001–4000 yr) to avoid initial model drift. All data are linearly detrended before analysis. Characteristics of the model’s Antarctic bottom water and its relationship with the Weddell Gyre have previously been described (e.g., Zhang and Delworth 2016). The impact of multidecadal Atlantic meridional overturning circulation (AMOC) variations on the SO using this model was also described in Zhang et al. (2017b). The realism of the model characteristics as described in these previous studies provides some level of confidence that this model is an appropriate tool for studies of the model predictability of SO.

### b. Methods

We first use the potential predictability variance fraction (Boer 2004; Boer and Lambert 2008) to give general information about the high predictability regions. Boer (2004) suggested that the total climate variability can be decomposed into the slow time scale “potentially predictable” component and unpredictable climate noise . The ppvf is therefore defined as a fraction of long time scale (or low frequency) variability with respect to the total variability . The term is the variance of *m*-yr mean SST, where *m* can be selected as any integer number. The high ppvf regions identify those areas in which long time scale variability stands out clearly from short time scale variability, and thus variability in these regions may be at least potentially predictable.

We then use the APT method to derive leading predictable components over these high ppvf regions. A standard measure of predictability (DelSole and Tippett 2009b) is defined as

where is climatological variance and is the ensemble forecast variance at lead time . This measure is close to one for a perfect forecast and close to zero when the ensemble forecast spread approaches the climatological spread.

The APT is defined as the integral of predictability over all lead times:

It is an integral measure of predictability and thus is independent of lead time. To maximize APT, we seek an inner product **q**^{T}**x**, where **q** is a projection vector, **x** is the state vector, and superscript T denotes the transpose operation. The component **q**^{T}**x** has respective forecast and climatological variances:

DelSole and Tippett (2009b) and Jia and DelSole (2011) pointed out that maximizing APT leads to an eigenvalue problem:

Since the control run data we used here only have a single ensemble member, a linear regression model is adopted to estimate APT. The regression model is written as

where denotes the predictor at time , is the predictand at time , is the regression coefficient at time, and is the residual term. The climatological and forecast matrices thus have the form of

where is the time-lagged covariance matrix and is the climatological variance. Substituting (7) into the eigenvalue problem (5) gives

The left term in (8) represents the integration of signal covariance, while the right term represents the total climatological covariance in which and denote the eigenvalue and projection vector, respectively.

When we apply this method to our control simulation, both the predictors and predictands are projected on the leading 30 principal components (PCs). The resulting 3000-yr length PCs are then split in half, as also seen in Jia and DelSole (2011). The first 1500 years of data from the control run, called training data, are used to maximize APT in Eq. (8), and the second 1500 years are kept for verification. As suggested by DelSole and Tippett (2009b), we use the squared multiple correlation to estimate the potential predictability; can represent the amount of variation in the predictand that is accounted for by the variation in the predictors and has a form of

The term **q** is calculated from the training data, while the covariance terms and are obtained from verification data. In general, the slower the decrease of with lead time, the larger the potential predictability. The statistical significance of APT is examined by Monte Carlo experiments. We generate two independent random matrices that have zero mean and unit variance and apply them to Eq. (8) to produce an ordered sequence of optimized APT values. This procedure is then repeated 100 times and the 95% eigenvalue from (8) was selected. If the APT value computed from training data exceeds the 95% value from the Monte Carlo experiments, the null hypothesis (white noise, unpredictable) will be rejected and the APT value from training data is significant at a 5% significance level.

## 3. Potential predictability over the Southern Ocean

### a. High predictability regions

Figure 1 shows the ppvf of 5-yr, 11-yr, and 25-yr mean SST over the global oceans in the CM2.1 model. In agreement with previous studies (e.g., Boer 2004; Boer and Lambert 2008), high-latitude regions exhibit relatively higher potential predictability than the middle and low latitudes. This contrast becomes more obvious when the averaging scale increases from 5 to 11 years (cf. Figs. 1a and 1b). The potential predictability of 11-yr mean SST is primarily concentrated in the North Atlantic, North Pacific, and SO, with comparable magnitudes in both hemispheres (Fig. 1b). When we consider the 25-yr mean SST, the SO potential predictability is even higher than the North Atlantic and North Pacific Oceans (Fig. 1c). The large values of ppvf over the SO in GFDL model indicate that the long time scale variability over the SO is very pronounced.

### b. APT analysis of SO SST

We identify the leading predictable components of SO SST using standard APT analysis. The analyzed area consists of all ocean grid points south of 35°S. Note that the results are not sensitive to the northern boundary choices, as long as the latitude is within the Southern Hemisphere (not shown). The leading two components have significant APT values at 5% significance level. The SST spatial pattern associated with the most predictable component (APT1) has loadings of the same sign over the SO, with maximum anomalies over the Amundsen–Bellingshausen–Weddell Seas (Fig. 2a). The APT value for this mode is 20.6 years. The corresponding time series of APT1 shows prominent multidecadal fluctuations (Fig. 2b), with a broad spectral peak around 70–120 years (Fig. 2c). The squared multiple correlation *R*^{2} of the leading predictable component as a function of time lag derived from independent control run is further shown in Fig. 2d. The *R*^{2} above the 95% confidence level denotes significant predictability. This figure shows that the APT1 mode has potential predictability up to 20 years. The traditional and damped persistence forecasts, which assume that the forecast equals the initial condition and the forecast decays exponentially in time, respectively, are shown in Fig. 2d as well. It can be seen that the skill arising from APT maximization is higher than both persistence forecasts.

The second most predictable component (APT2) of SO SST is characterized by a dipole structure, with SST anomalies of one sign over the Weddell Sea and SST anomalies of the opposite sign in the Amundsen–Bellingshausen Seas (Fig. 3a). The associated time series has a pronounced multidecadal oscillation, but is quite noisy compared to APT1 (Fig. 3b vs Fig. 2b). The power spectrum of APT2 time series reveals a quasi-70–120-yr peak (Fig. 3c) that also appears in APT1 (Fig. 2c). Figure 3c also shows substantial variances are distributed in the 2–50-yr frequency bands in the APT2 power spectrum, in sharp contrast to APT1 (Fig. 2c). This leads to a relatively small variance fraction in 70–120-yr frequency bands and thus a noisy APT2 time series. The *R*^{2} of this second most predictable component indicates a potential predictability up to 5 years (Fig. 3d), which is much shorter than the first predictable mode due to noisy characteristics. The APT2 predictability is only slightly higher than the persistence forecasts (Fig. 3d), suggesting that the skill mainly arises from the SST persistence.

The coherent spectrum of APT1 and APT2 time series shows high coherences over their common peak period 70–120 yr (Fig. 4a), suggesting that the leading two predictable components may have the same ocean origin. To confirm our hypothesis, we conduct a lead–lag correlation analysis between these two time series (Fig. 4b). As expected, the simultaneous correlation is zero due to the orthogonality of the APT decomposition. Significant positive (negative) correlations are found when the APT1 leads (lags) APT2 by 10–30 years. These lead and lag times account for approximately a quarter of the APT1/APT2 period (70–120 yr). These analyses imply that the two predictable components appear to vary in quadrature.

## 4. Ocean origin of high decadal predictability of SO SST

### a. Climate fluctuations associated with leading predictable modes

To understand the physical processes associated with the leading predictable components, we regress several important variables onto the APT1 and APT2 time series, respectively (Figs. 5 and 6). Figure 5a exhibits the surface net heat flux and mixed layer depth (MLD) anomalies associated with the APT1 time series. The SO experiences broad negative heat flux anomalies that tend to damp positive surface temperature anomalies. This implies that the uniform SO SST warming in APT1 originates from the ocean dynamics, instead of atmosphere forcing. The MLD change shown in Fig. 5a displays a strong positive anomaly over the Weddell Sea, indicating strong deep convection there. Note that the long-term mean global meridional overturning circulation (GMOC) has a negative value south of 60°S, which represents an anticlockwise cell that denotes the strength of Antarctic Bottom water (AABW) formation as well as deep convection (Fig. 7a). Figure 5b shows prominent negative GMOC anomalies south of 20°S, suggesting a strengthening and northward extension of the AABW cell. In the mean state the subsurface is warmer than the surface in the region of the SO. Therefore, the spinup of AABW cell drives a subsurface–surface temperature dipole in the SO, with a cooling anomaly in the subsurface and a warming anomaly at the surface (Fig. 5c) that corresponds to a decrease of Antarctic sea ice (Fig. 5d). The surface wind is characterized by a zonally oriented anticyclone around 40°–60°S band, which is very likely due to local SST feedback (Zhang et al. 2017b). The easterly wind anomaly around 40°S also favors warm SST in the midlatitude due to anomalous warm southward Ekman transport. The heat flux and MLD anomalies associated with the APT2 time series show opposite signs with the APT1 (cf. Figs. 6a and 5a), suggesting a weakening of deep convection over the Weddell Sea. Accordingly, the GMOC anomaly shows a spindown of AABW cell (Fig. 6b). Compared to APT1, the GMOC change is relatively weak and mainly confined south of 60°S (cf. Figs. 5b and 6b). The associated zonal mean temperature shows a weak cold surface–warm subsurface dipole structure over the SO (Fig. 6c). In contrast to the uniform sea ice response in APT1, the sea ice change associated with APT2 exhibit a dipole pattern, with sea ice increase in the Weddell Sea and sea ice decrease over the Amundsen–Bellingshausen Seas (Fig. 6d). The sea ice anomaly over the Amundsen–Bellingshausen Seas is not only related to the SST anomalies but also linked with the surface wind. As shown in Fig. 6d, there is an anticyclonic wind around 160°–40°W over the SO, which corresponds to a northwest wind anomaly over the Amundsen–Bellingshausen Seas. The northwest wind favors poleward warm temperature advection and thus a decrease of sea ice there.

The above regression analyses suggest that the leading two predictable components of SO SST are very likely to be associated with deep convection changes. To test this hypothesis, we examine the SO deep convection characteristics in the CM2.1 model. As mentioned above, we use the AABW cell anomaly to represent the SO deep convection fluctuations. The strength of the AABW cell each year is defined as the minimum value of the streamfunction south of 60°S (Fig. 7a). Note that a negative AABW cell anomaly indicates a stronger than usual overturning cell. The time series of the AABW cell index (Fig. 7b) has pronounced multidecadal variability at time scales of approximately 70 to 120 years (Fig. 7c); this coincides with the broad spectral peaks of the APT time series (Fig. 7c vs Figs. 2c and 3c). We also show in Fig. 7d the lead lag correlation between the AABW cell index and the APT1 time series. It shows a negative correlation as low as −0.6 when the AABW leads the APT1 by about 5 years. Since the APT1 and APT2 time series are in quadrature, significant correlations are also found between the AABW index and APT2 with some time lags (not shown).

### b. Southern Ocean multidecadal variability

To further confirm the close relationship between the leading predictable components and SO deep convection, we show in Figs. 8a–e the multidecadal cycle of AABW cell. The AABW cell cycle is obtained by the lagged regression of GMOC anomalies upon the AABW cell index. At a lag of 0 yr, the AABW cell is in its mature positive phase, with a maximum increase south of 60°S and a northward extension to 40°S (Fig. 8a). As we move forward from lag 0, the GMOC anomalies south of 60°S gradually weaken, while the northward extension becomes stronger and stronger (Fig. 8b). At a lag of 20 yr, a positive GMOC anomaly emerges south of 60°S. This positive GMOC anomaly then intensifies and gradually spreads northward, which in turn weakens the negative GMOC anomaly in the north (Figs. 8c,d). At a lag of 40 yr, the AABW phase is totally flipped and reaches its mature negative phase (Fig. 8e). A close examination finds that the spatial structure of quasi mature phase of AABW cell (Fig. 8a) closely resembles the GMOC anomalies associated with the APT1 (Fig. 5b). Similarly, the transition phase of AABW cycle (Fig. 8c) matches with the GMOC anomalies associated with the APT2 very well (Fig. 6b).

We show in Figs. 8f–j the multidecadal SST cycle associated with the deep convection. During the AABW cell mature positive phase, the SO experiences broad warming anomalies, with maximum values over the Weddell Sea (Fig. 8f). The warm SST over the SO corresponds to a zonally oriented anticyclone wind, with easterly anomalies around 40°S and westerly anomalies at 75°S. We note that the mature phase of the SO SST cycle here is in good agreement with the SST pattern in APT1 (cf. Figs. 8f and 2a). Accompanied with the AABW cell weakening south of 60°S (Figs. 8a,b), the positive SST anomalies over the Weddell Sea gradually weaken (Figs. 8f,g). At the same time, the southeast Pacific SST warming gradually spreads to the equator through the fast positive wind–evaporation–SST (WES) feedback and slow weakening of subtropical cell (e.g., Ma and Wu 2011) (Figs. 8f–h). The warm SST anomaly over the tropical Pacific further induces a positive phase of the Pacific–South American (PSA) teleconnection (Mo and Higgins 1998). The PSA teleconnection corresponds to a wavenumber 3 in the midlatitudes with a high geopotential height and an anticyclonic circulation over the Amundsen–Bellingshausen Seas that favors warm poleward advection and thus warm SST there (Fig. 8h). At a lag of 20 yr, a cooling SST anomaly appears in the Weddell Sea (Fig. 8h) resulting from the emergence of a weak AABW cell shown in Fig. 8c. The SST over the SO is characterized by a dipole pattern at this moment, with a warm SST in the Amundsen–Bellingshausen Seas and a cold SST in the Weddell Sea. The negative SST anomaly in the Weddell Sea further grows in the Weddell Sea and then extends to the entire SO (Figs. 8h-j). Apparently, the SST anomalies over the Amundsen–Bellingshausen Seas lag the SST anomalies over the Weddell Sea due to slow advection times. At lags of 40 yr, the SO is almost covered by the negative SST anomalies (Fig. 8j), which reach to the opposite phase of deep convection. We note again that the transition phase of SO SST cycle matches very well with the SST pattern in APT2 (cf. Figs. 8h and 3a). These SST pattern similarities suggest that the leading two predictable components of SO SST originate from the internal multidecadal cycle of SO deep convection. The first component arises from the quasi-mature phase of deep convection, while the second component is contributed from the transition phase of deep convection.

The associated sea ice and subsurface temperature variabilities (Fig. 9) are physically consistent with our previous analyses. The sea ice primarily follows the SST changes, with a cold (warm) SST anomaly corresponding to a sea ice increase (decrease). Thus, the mature positive phase sea ice at lag 0 yr is characterized by a sea ice decrease over the SO (Fig. 9a), whereas the transition phase sea ice at lag 25 yr exhibits a sea ice decrease in the Amundsen–Bellingshausen Seas and a sea ice increase in the Weddell Sea (Fig. 9c). Figures 9f–j show the multidecadal zonal mean temperature cycle. As expected, the temperature response agrees with the deep convection change (cf. Figs. 9f–j and 8a–e). The spinup (spindown) of AABW cell brings subsurface warm water to the surface, thereby leading to a warm (cold) SST in the surface and a cold (warm) temperature in the subsurface. The dipole temperature weakens as the AABW cell spins down and vice versa. Again, these sea ice and zonal mean temperature anomalies in the mature and transition phases match with the responses in APT1 and APT2, respectively (cf. Figs. 9a,f and 5c,d; cf. Figs. 9c,h and 6c,d).

We note that the most predictable SST (APT1) time series over the SO lags the AABW cell index by about 5 years (Fig. 7d). The delayed SST response primarily arises from the slow adjustment of the ocean that consists of advection/wave propagation (Zhang and Delworth 2016), which is also seen in the CM2.1 fully coupled control run. Figure 10a exhibits the SO (50°–75°S, 0°–360°E) area averaged SST time series and the Weddell Sea (75°–55°S, 52°W–0°E) area averaged SST time series as well as the AABW cell index. Their lead–lag correlations are shown in Fig. 10b. All three indices have pronounced multidecadal fluctuations and they are highly correlated. The AABW cell index is simultaneously correlated with the local (Weddell Sea) SST due to strong deep convection there. In contrast, the maximum correlation between the AABW cell index and remote SO SST occurs when the AABW cell leads by about 5 years (Fig. 10b).

### c. Mechanisms contributing to SO multidecadal variability

The periodic strengthening and weakening of deep convection caused by subsurface advection heating and surface freshening results in model multidecadal oscillation (Figs. 11 and 12). These mechanisms have much in common with similar variability found in the Kiel Climate Model (Martin et al. 2013). Figure 11a shows the time evolution of annual mean vertical temperature anomaly averaged over the Weddell Sea. During active convection (strong AABW cell anomaly in Fig. 11c), the temperature distribution is almost homogeneous over the entire water column, corresponding to weak column stability, *N*^{2} (buoyancy frequency, ; Fig. 11d). During weak convection (weak AABW cell anomaly in Fig. 11c), the ocean becomes static stable with large values of *N*^{2}. At this time, the heat tends to accumulate at middepth, while the SST cools because the surface water directly contacts with the atmosphere and is strongly damped by the surface heat flux The middepth heat spreads over time, warms the entire water column below 300 m (Fig. 11a), destabilizes the ocean stratification from below (Fig. 11e), leads to a decrease of *N*^{2} (Fig. 11d), and eventually triggers the occurrence of deep convection. The convection weakening is preceded by subsurface heat depletion as well as surface freshening (Figs. 11a,b). The stabilizing freshwater cap gets thicker and fresher as convection weakens. The minimum surface salinity (Fig. 11b), and therefore the strongest stabilizing effect (Fig. 11d), occurs just slightly before the peak phase of weak convection. Obviously, the increase of static stability *N*^{2} from strong convection to weak convection comes from the salinity contribution (Figs. 11d–f).

We now focus on what are the main processes responsible for subsurface warming and surface freshening. Figure 12a shows the temperature budget at 2000 m for the same period and region shown in Fig. 11a. It can be seen that the buildup of heat in the subsurface is mainly due to temperature advection. The diffusive warming also contributes positively, but with a much smaller magnitude. Further decomposition finds that the advective warming is dominated by horizontal temperature advection, while the vertical advection contributes negligibly (Fig. 12b). The horizontal temperature advection is largely associated with the horizontal Weddell Gyre. Before the weakening of deep convection, the westward return flow in the southern branch of the Weddell Gyre effectively drags midlatitude heat (mainly the North Atlantic Deep Water) into the Weddell Sea (Fig. 12c). This strong barotropic clockwise gyre exists over almost the entire water column and its strength is strongly associated with the deep convection itself due to interactions between the AABW outflow and topography (Zhang and Delworth 2016). The gyre still exists when the convection spins down, albeit with a weak amplitude (Fig. 12d). This guarantees a persistent inverted temperature structure over the SO, which is necessary but not sufficient for an oscillation.

Figures 12e and 12f show the surface salinity budget over the Weddell Sea. Apparently, the surface freshening before the onset of weak convection is caused by salinity advection and surface freshwater forcing as well as diffusion (Fig. 12e). The contribution from diffusion is relatively small compared to the other two factors. Figure 12f further shows the salinity advection mainly arises from horizontal advection. The Weddell Sea surface is salty during active convection because of high-salinity water from the subsurface. The salinity over adjacent basins is relatively low, thus the horizontal salinity advection from remote regions can induce freshening over the Weddell Sea. The positive contribution of surface salt flux mainly comes from the sea ice melting in response to warm SST during active convection (Fig. 12f).

In brief, the SO multidecadal variability mainly a result of halocline and inverted thermocline structure over the subpolar region. Horizontal warm advection at depth due to the Weddell Gyre warms up the subsurface water and induces convection (red color in Fig. 13). Anomalous freshwater advected from remote regions and sea ice melting act as a fresh cap and thus weaken convection (blue color in Fig. 13). The whole cycle repeats itself. The time scale is determined by the rate of subsurface heating and surface freshening. In the CM2.1 model, the recharge and discharge processes of subsurface heat reservoir are related to ocean advection, which are quite slow and thus form low-frequency oscillation. Although the atmosphere response can somehow change the amplitude of SST and other related variables (e.g., Figs. 8f–j), it is not the main driver and cannot determine the long time period.

## 5. Climate impacts

In this section, we examine the potential multiyear predictability of surface air temperature (SAT) and precipitation over the Antarctic continent. We find that the multiyear predictability of land variables using the land predictor itself is lower than the predictability using global SST (not shown). This suggests that the land predictability on interannual-to-decadal time scales is primarily driven by SST (Hoerling and Kumar 2003; Held et al. 2005). Thus, we use a generalized APT method (GAPT) (Jia and DelSole 2011), which is similar to the standard APT described in section 2, except that the predictor and predictand are two different variables. Here the predictor is global SST, while the predictand is SAT or precipitation.

We show in Fig. 14 the most predictable component of SAT over the Antarctic continent. The SAT physical pattern has a uniform warming over the entire Antarctic land area, with maximum amplitudes over the Antarctic Peninsula (Fig. 14a). The *R*^{2} values in the verification data suggest that the Antarctic SAT is able to be predicted 6 years in advance (Fig. 14c). The SAT time series has a pronounced 70–110-yr peak (Fig. 14b), implying a potential linkage with the SO deep convection. The SST regression pattern associated with the APT1 time series shows notable SST warming over the SO, with negligible signals over the Northern Hemisphere (Fig. 14d). The maximum SST warming occurs over the Weddell Sea where deep convection dominates. This SST pattern highly resembles the mature phase SST associated with the deep convection fluctuations (cf. Figs. 14d and 9a). Thus, we conclude that the most predictable SAT over the Antarctic continent results from the SO SST memory that is controlled by the deep convection activity.

Similar to the SAT, the most predictable component of precipitation over the Antarctic continent primarily arises from the prediction skill of SO SST that is closely linked with SO deep convection. The physical pattern of precipitation displays positive anomalies over land where the adjacent SO has significant SST warming (cf. Figs. 15a and 15d). The power spectrum of the GAPT1 time series, again, shows a similar frequency peak with the SO SST and AABW cell (Fig. 15b vs Figs. 2c and 7d). The *R*^{2} values of precipitation are lower than that of SAT due to the noisy characteristics of precipitation (cf. Figs. 14c and 15c). However, the potential predictability of precipitation can still be up to 4 years (Fig. 15c).

## 6. Summary, discussion, and conclusions

By taking advantage of the GFDL CM2.1 4000-yr control run integration, we investigate the potential decadal predictability of SO SST in the present paper. We use a new statistical optimization technique, called APT analysis (DelSole and Tippett 2009a,b), to identify the leading predictable components of SO SST on decadal time scales. The APT analysis maximizes an integrated prediction variance obtained from a linear regression model, in which both the predictor and predictand are SST. Note that the long-term control integration does not include anthropogenic forcing or changes in natural forcing from volcanoes or interannual variations of solar irradiance. The potential predictability shown here is therefore purely from internal variability.

The most predictable component of SO SST can be predicted in an independent verification data by a linear regression model, with significant skill up to 20 years. The predictable pattern has a uniform SST sign over the SO, with maximum values over the Amundsen–Bellingshausen–Weddell Seas. The associated APT time series has a 70–120-yr spectral peak. This predictable pattern is closely related to the mature phase of the SO internal variability that originates from deep convection fluctuations. In the CM2.1 model, deep convection mainly occurs over the Weddell Sea and has multidecadal fluctuations on a 70–120-yr time scale. This multidecadal time scale selection is largely associated with the recharge processes of heat reservoir in the deep ocean. Slow subsurface ocean processes provide long time scales that give rise to decadal predictability of SO SST. The SO SST has significant climate impacts on the SAT and precipitation over the Antarctic continent. The SAT and precipitation can be potentially predictable up to 6 yr and 4 yr in advance, respectively. These multiyear prediction skills arise from the SO SST, which is again attributed to the internal deep convection fluctuations over the Weddell Sea.

The second most predictable component of SO SST is characterized by a dipole structure, with SST anomalies of one sign over the Weddell Sea and SST anomalies of the opposite sign over the Amundsen–Bellingshausen Seas. This component has statistically significant prediction skill for 6 years based on a linear regression model. A close examination reveals that this dipole mode primarily arises from the transition phase of the dominant pattern of SO internal variability. Again, the slow ocean memory associated with the SO deep convection provides the multiyear prediction skill of this second most predictable component. Interestingly, this second component corresponds to a sea ice dipole structure over the Amundsen–Bellingshausen–Weddell Seas, which somewhat resembles the observed sea ice trend in recent years (e.g., Li et al. 2014). The associated surface wind over the SO characterized by a cyclone (or anticyclone) around 160°–40°W favors sea ice dipole formation, which is also consistent with what was found in observation. These similarities provide a hypothesis that some prominent trends observed during the recent decades in the Southern Hemisphere may have some contributions from internal variability in the SO that is strongly associated with deep convection fluctuations.

To provide some perspective on the robustness of our results, we perform the same APT diagnostics on a long control integration of a different GFDL climate model, CM3 (Donner et al. 2011). The CM3 model has an ocean component that is quite similar to CM2.1, but the atmospheric component of CM3 has substantial differences from CM2.1. Similar to CM2.1, the most predictable SST pattern in CM3 displays a uniform sign over the SO, while the second most predictable SST pattern shows a dipole structure (Figs. 16a,b). The maximum SST centers associated with these two modes are primarily over the continental shelf regions of the Weddell and Ross Seas (Figs. 16a,b), which are somewhat different from CM2.1 (Figs. 2 and 3). These SST center differences are largely associated with the deep convection position in the two models (not shown). The SO deep convection mainly occurs in the Weddell Sea in the CM2.1 model, including both over the open ocean and the continental shelf. In contrast, the deep convection in CM3 model takes places in continental shelf regions of the Weddell and Ross Seas.

The power spectrum of APT1 and APT2 time series in CM3 model shows prominent spectral peaks around 300 years, which are longer than that in CM2.1 model (Fig. 16c vs Figs. 2c and 3c). This leads to a long persistence time of SST and therefore a long predictability skill, as presented in Figs. 16d and 16e. The predictive skill can be up to 30 years for APT1 component and 10 years for APT2 mode. In agreement with that in the CM2.1 model, these leading predictable components are found to be closely linked with the SO deep convection fluctuations (Fig. 16f). The period difference in these two models seems to arise from differences in the warming rate in the subsurface (not shown) that are related to the Weddell Gyre strength and warming magnitude at depth in the midlatitude. We will examine this issue further in future work.

Our diagnostic approach for the decadal predictability of SO SST in the current paper suggests that if we could correctly initialize the SO deep convection in the numerical forecast model, the future evolution of SO SST and its associated climate impacts might be predictable on decadal scales. Such predictions would ideally be performed using models with simulations of the SO that are as realistic as possible. In addition, enhanced ocean observations, particularly subsurface observations over the far SO, are also needed to characterize the state of the SO. An important caveat is that the realism of model’s simulation of the SO will impact how relevant such potential predictive skill is for predictions of the real climate system. The decadal prediction skill of SO SST based on real decadal hindcasts/forecasts is currently under investigation and will be the topic of a forthcoming paper (Zhang et al. 2017a).

## Acknowledgments

We are grateful to Soyna Legg and Carolina O. Dufour for their constructive comments on an early version of the paper, as well as three anonymous reviewers who provided extremely insightful and valuable feedback and suggestions. The work of T. Delworth is supported as a base activity of NOAA’s Geophysical Fluid Dynamics Laboratory. L. Zhang and L. Jia are supported through Princeton University under block funding from NOAA/GFDL.

## REFERENCES

## Footnotes

© 2017 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).