Interannual drivers of the seasonal cycle of CO 2 in the Southern Ocean

. Resolving and understanding the drivers of variability of CO 2 in the Southern Ocean and its potential climate feedback is one of the major scientiﬁc challenges of the ocean-climate community. Here we use a regional approach on empirical estimates of p CO 2 to understand the role that seasonal variability has in long-term CO 2 changes in the Southern Ocean. Machine learning has become the preferred empirical modelling tool to interpolate time- and location-restricted ship measurements of p CO 2 . In this study we use an ensemble of three machine-learning products: support vector regression (SVR) and random forest regression (RFR) from Gregor et al. (2017), and the self-organising-map feed-forward neural network (SOM-FFN) method from Landschützer et al. (2016). The interpolated estimates of (cid:49)p CO 2 are separated into nine regions in the Southern Ocean deﬁned by basin (Indian, Paciﬁc, and Atlantic) and biomes (as deﬁned by Fay and McKinley, 2014a). The regional approach shows that, while there is good agreement in the overall trend of the products, there are periods and regions where the conﬁdence in estimated (cid:49)p CO 2 is low due to disagreement between the products. The regional breakdown of the data highlighted the seasonal decoupling of the modes for summer and winter interannual variability. Winter interannual variability had a longer mode of variability compared to summer, which varied on a 4–6-year timescale. We separate the analysis of the (cid:49)p CO 2 and its drivers into summer and winter. We ﬁnd that understanding the variability of (cid:49)p CO 2 and its drivers on shorter timescales is critical to resolving the long-term variability of (cid:49)p CO 2 . Results show that (cid:49)p CO 2 is rarely driven by thermodynamics during winter, but rather by mixing and stratiﬁcation due to the stronger correlation of (cid:49)p CO 2 variability with mixed layer depth. Summer p CO 2 variability is consistent with chlorophyll a variability, where higher concentrations of chlorophyll a correspond with lower p CO 2 concentrations. In regions of low chlorophyll a concentrations, wind stress and sea surface temperature emerged as stronger drivers of (cid:49)p CO 2 . In summary we propose that sub-decadal variability is explained by summer drivers, while winter variability contributes to the long-term changes associated with the SAM. This approach is a useful framework to assess the drivers of (cid:49)p CO 2 but would greatly beneﬁt from improved estimates of (cid:49)p CO 2 and a longer time series.

Abstract. Resolving and understanding the drivers of variability of CO 2 in the Southern Ocean and its potential climate feedback is one of the major scientific challenges of the ocean-climate community. Here we use a regional approach on empirical estimates of pCO 2 to understand the role that seasonal variability has in long-term CO 2 changes in the Southern Ocean. Machine learning has become the preferred empirical modelling tool to interpolate time-and locationrestricted ship measurements of pCO 2 . In this study we use an ensemble of three machine-learning products: support vector regression (SVR) and random forest regression (RFR) from , and the self-organising-map feedforward neural network (SOM-FFN) method from Landschützer et al. (2016). The interpolated estimates of pCO 2 are separated into nine regions in the Southern Ocean defined by basin (Indian, Pacific, and Atlantic) and biomes (as defined by Fay and McKinley, 2014a). The regional approach shows that, while there is good agreement in the overall trend of the products, there are periods and regions where the confidence in estimated pCO 2 is low due to disagreement between the products. The regional breakdown of the data highlighted the seasonal decoupling of the modes for summer and winter interannual variability. Winter interannual variability had a longer mode of variability compared to summer, which varied on a 4-6-year timescale. We separate the analysis of the pCO 2 and its drivers into summer and winter. We find that understanding the variability of pCO 2 and its drivers on shorter timescales is critical to resolving the long-term variability of pCO 2 . Results show that pCO 2 is rarely driven by thermodynamics during winter, but rather by mixing and stratification due to the stronger correlation of pCO 2 variability with mixed layer depth. Summer pCO 2 variability is consistent with chlorophyll a variability, where higher concentrations of chlorophyll a correspond with lower pCO 2 concentrations. In regions of low chlorophyll a concentrations, wind stress and sea surface temperature emerged as stronger drivers of pCO 2 . In summary we propose that sub-decadal variability is explained by summer drivers, while winter variability contributes to the long-term changes associated with the SAM. This approach is a useful framework to assess the drivers of pCO 2 but would greatly benefit from improved estimates of pCO 2 and a longer time series.

Introduction
The Southern Ocean plays a key role in the uptake of anthropogenic CO 2 (Khatiwala et al., 2013;DeVries et al., 2017). Moreover, it has been shown that the Southern Ocean is sensitive to anthropogenically influenced climate variability, such as the intensification of the westerlies (Le Quéré et al., 2007;Lenton et al., 2009;Swart and Fyfe, 2012;DeVries et al., 2017). Until recently, the research community has not been able to quantify the contemporary changes of CO 2 in the Southern Ocean accurately due to a paucity of observations, let alone understand the drivers (Bakker et al., 2016). Empirical models provide an interim solution to this challenge until prognostic ocean biogeochemical models are able to represent Southern Ocean CO 2 fluxes adequately Rödenbeck et al., 2015;Mongwe et al., 2016).
The research community agrees on large changes in CO 2 fluxes in the Southern Ocean from a weakening sink in the 1990s to a strengthening sink in the 2000s; however, there is disagreement over the drivers of the changes in CO 2 uptake (Lovenduski et al., 2008;Landschützer et al., 2015;DeVries et al., 2017;Ritter et al., 2017). This study aims to understand the drivers of the changing CO 2 sink in the Southern Ocean, based on an ensemble of empirical estimates using a seasonal analysis framework.
Empirical methods estimate CO 2 by extrapolating sparse ship-based CO 2 measurements using proxy variables. The proxies are often observable by satellite but may include climatologies or output from assimilative models. Empirical methods have improved our understanding of CO 2 trends in the Southern Ocean by increasing the data coverage. However, there is still disagreement between many of the methods due to the paucity of data and the way in which each method interpolates sparse data Ritter et al., 2017).
In a key study, Landschützer et al. (2015) showed, using an artificial neural network (ANN), that there was significant strengthening of Southern Ocean CO 2 uptake during the period 2002-2010. While previous studies suggested that changes in wind strength have led to changes in meridional overturning and thus CO 2 uptake (Lenton and Matear, 2007;Lovenduski et al., 2007;Lenton et al., 2009;DeVries et al., 2017), Landschützer et al. (2015) suggested that atmospheric circulation has become more zonally asymmetric since the mid-2000s, which has led to an oceanic dipole of cooling and warming. The net impact of cooling and warming, together with changes in the DIC and TA (dissolved inorganic carbon and total alkalinity), led to an increase in the uptake of CO 2 . During this period, southward advection in the Atlantic basin reduced upwelled DIC in surface waters, overcoming the effect of the concomitant warming in the region. Conversely, in the eastern Pacific sector of the Southern Ocean, strong cooling overwhelmed increased upwelling . This is supported by observations from the Drake Passage and south of Australia showing that variability of upwelling has affected pCO 2 (Munro et al., 2015;Xue et al., 2015).
In a subsequent study, Landschützer et al. (2016) proposed that interannual variability of CO 2 in the Southern Ocean is tied to the decadal variability of the Southern Annular Mode (SAM) -the dominant mode of atmospheric variability in the Southern Hemisphere (Marshall, 2003). This concurs with previous studies, which suggested that the increase in the SAM during the 1990s resulted in the weakening of the Southern Ocean sink (Le Quéré et al., 2007;Lenton and Matear, 2007;Lovenduski et al., 2007;Lenton et al., 2009;Xue et al., 2015). The work by Fogt et al. (2012) bridges the gap between the proposed asymmetric atmospheric circulation of Landschützer et al. (2015) and the observed correlation with the SAM of Landschützer et al. (2016). Fogt et al. (2012) show that changes in the SAM have been zon-ally asymmetric and that this variability is highly seasonal, thus amplifying or suppressing the amplitude of the seasonal mode.
Assessing the changes through a seasonal framework may thus help shed light on the drivers of CO 2 in the Southern Ocean. Southern Ocean seasonal dynamics suggest that the processes driving pCO 2 are complex, but with two clear contrasting extremes. In winter, the dominant deep mixing and entrainment processes are zonally uniform, driving an increase in pCO 2 , with the region south of the Polar Front (PF) becoming a net source and weakening the net sink north of the PF . In summer, the picture is more spatially heterogeneous, with net primary production being the primary driver of variability (Mahadevan et al., 2011;Thomalla et al., 2011;Lenton et al., 2013). The competing influence between light and iron limitation results in heterogeneous distribution of chlorophyll a (Chl a) in both space and time, with similar implications for pCO 2 ( Thomalla et al., 2011;Carranza and Gille, 2015). The interaction between the large-scale drivers, such as wind stress, surface heating, and mesoscale ocean dynamics, are the primary cause of this complex picture (McGillicuddy, 2016;Mahadevan et al., 2012). Some regions of elevated mesoscale and sub-mesoscale dynamics, mainly in the Sub-Antarctic Zone (SAZ), are also characterized by strong intra-seasonal modes in summer primary production and pCO 2 (Thomalla et al., 2011;Monteiro et al., 2015). In general, the opposing effects of mixing and primary production result in the seasonal cycle being the dominant mode of variability in the Southern Ocean .
In this study we examine winter and summer interannual variability of pCO 2 in the Southern Ocean between 1998 and 2014 to understand the drivers of long-term changes in CO 2 uptake.

Empirical methods and data
In this study we use three machine-learning methods: random forest regression (RFR), support vector regression (SVR) and a self-organising-map feed-forward neural network (SOM-FFN). RFR and SVR are introduced in  and SOM-FFN is presented in Landschützer et al. (2014). In brief, the RFR approach is an ensemble of decision trees that provides non-linear regression by combining many high variance-low bias estimators . SVRs are in principle similar to a single-hidden-layer FFN, except that SVR statistically determines the complexity of the problem, which is analogous to the hidden layer structure that is typically determined heuristically. The SOM-FFN method is a two-step neural network approach that first clusters data (SOM) and then applies a regression model (FFN) to each cluster.
The SVR and RFR implementations used in this study are trained with the monthly 1 by 1 • gridded SOCAT (Surface Ocean CO 2 Atlas) v3 dataset (Bakker et al., 2016). The SOM-FFN (run ID: netGO5) used in this study was trained with SOCAT v4 . Table 1 shows the proxy variables used for each of the methods. Sea surface salinity (SSS) and mixed layer depth (MLD) for SVR and RFR are from Estimating the Circulation and Climate of the Ocean, Phase II (ECCO 2 ) (Menemenlis et al., 2008). The use of these assimilative modelled products may in some cases produce results that are unrealistic. This may have influenced the use of the de Boyer Montégut et al. (2004) MLD climatology in the SOM-FFN, where ECCO 2 was used in previous iterations of the product. The trade-off of using the climatology is that no interannual changes in MLD are taken into account. We acknowledge that using different proxy variables could result in different pCO 2 estimates, but comparing the different proxies used in each of the CO 2 products is beyond the scope of this study. Other data sources that are consistent between methods are sea surface temperature (SST) and sea-ice fraction by Reynolds et al. (2007), Chl a by Maritorena and Siegel (2005), absolute dynamic topography (ADT) generated by DUACS and distributed by AVISO (ftp://ftp.aviso.oceanobs.com, last access: 12 July 2012), and xCO 2 (CDIAC, 2016) with pCO 2(atm) calculated from interpolated xCO 2 using NCEP2 sea level pressure (Kanamitsu et al., 2002). In the case of Chl a for SVR and RFR,  filled the cloud gaps with climatological Chl a. Note that ADT coverage is limited to regions of no to very low concentrations of sea-ice cover, so estimates for SVR and RFR do not extend into the ice-covered regions during winter. Our analyses are thus limited to the region north of the maximum sea-ice extent.
Seasonality of the data is preserved by transforming the day of the year (j ) and is included in both SVR and RFR analyses: Transformed coordinate vectors are passed to SVR only using n-vector transformations of latitude (λ) and longitude (µ) (Gade, 2010;Sasse et al., 2013), with N containing the following: Wind speed, while not used in the empirical methods, is used in the assessment of the drivers of CO 2 . We use CCMP v2, which is an observation-based product that combines remote sensing, ship, and weather buoy data (originally published by Atlas et al., 2011 andupdated by Wentz et al., 2015). Swart et al. (2015a) compared a number of wind reanalysis products with CCMP v1 (where CCMP was the benchmark). The authors found that many of the reanalysis products had spurious trends, particularly in the Southern Hemisphere where data are sparse. Our choice of CCMP, which is based on observations, aims to minimise the assumptions that are otherwise made by reanalysis products.

Uncertainties
The machine-learning approaches used in this study are by no means able to estimate pCO 2 with absolute certainty. To account for the uncertainty, we use the same approach as Landschützer et al. (2014) to calculate total errors for each of the methods: where e m(t) is the total error associated with a method (m); e meas is the error associated with SOCAT measurements, which is fixed at 5 µatm (Pfeil et al., 2013), where 1 µatm is 101 325 Pa; e grid is the 5 µatm error associated with gridding the data into monthly by 1 • bins . Lastly e map is the root mean squared error (RMSE) calculated for each method, as shown in Table 1 taken from . These errors are used to calculate the average "withinmethod" error as defined by Gurney et al. (2004): where e m(t) is the method-specific error as defined in Eq.
(3) and M is the number of methods (three in this case). For a measure of the difference between methods we use the "between-method" approach used in Gurney et al. (2004): where S m is the method estimate of pCO 2 and S is the mean of the methods. This is analogous to the standard deviation (for a known population size). We later use an adaptation of this metric as a threshold to determine the confidence around anomalies.

Regional coherence framework
Southern Ocean CO 2 is spatially heterogeneous, both zonally and meridionally . In order to understand this heterogeneity, we used the three southernmost biomes defined by Fay and McKinley (2014a), as done in Rödenbeck et al. (2015). From north to south these are the subtropical seasonally stratified (STSS), sub-polar seasonally stratified (SPSS), and seasonally ice-covered region (ICE). These three biomes are comparable to the SAZ, PFZ (Polar-Frontal Zone) and MIZ (Marginal Ice Zone) respectively and will be Table 1. Three empirical methods used in the ensemble. RFR and SVR are described in . SOM-FFN is from Landschützer et al. (2016). SST = sea surface temperature, MLD = mixed layer depth, SSS = sea surface salinity, ADT = absolute dynamic topography, Chl a = chlorophyll a, pCO 2(atm) = fugacity of atmospheric CO 2 , xCO 2(atm) = mole fraction of atmospheric CO 2 , (lat., long.) = N−vector transformations of latitude and longitude, t(day of year) = trigonometric transformation of the day of the year. Note that SOM-FFN uses the de Boyer Montégut et al. (2004) climatology for MLD (dBM2004). The root mean squared errors (RMSEs) listed in the last column are for the Southern Ocean from .

Method
Input

Results and discussion
The first section of the results examines the uncertainties of the ensemble and its members. We then look at the seasonal cycle of the ensemble mean in time and space. This is done to lay the foundation for the interpretation of the results when assessed with the regional framework. In the regional interpretation the estimates are decomposed into nine regions, as shown in Fig. 1. Lastly, we implement a seasonal decomposition of the estimates to interpret the drivers of the changes observed in pCO 2 .

Ensemble member performance and variability
We use the RMSE scores as presented in  with abbreviated results shown in Table 1. The SOM-FFN method has the best score (14.84 µatm). SVR scores the lowest (24.04 µatm), but is still included due to the method's sensitivity to sparse data, which is favourable to the poorly sampled winter period . This compliments the RFR method, which scores well (16.45 µatm) but is prone to being insensitive to sparse data . These RMSE scores are used to calculate the total errors for each method and region using Eq. (3), where the measurement and mapping errors are both 5 µatm each (Pfeil et al., 2013;. These results are shown in Table 2. Total errors are used to calculate the within-method error, which is an estimate of the combined total error of the three machine-learning methods (Eq. 4). The between-method errors are the mean of the standard deviation between the methods (Eq. 5). The within-method errors are much larger than the between-method errors (Table 2). However, the withinmethod errors are normally distributed and are mechanistically consistent . The between-method The three biomes used in this study, SAZ, PFZ, and MIZ, are defined by Fay and McKinley (2014a). The regions are also split by basin. Fig. 2c) is thus used to determine whether observed variability is consistent between the three methods. Figure 2 shows the pCO 2 time series for the SAZ and PFZ. Note that we exclude the MIZ from the remaining analyses due to large E b and E w ( Table 2) and inconsistent coverage between products due to sea-ice cover (MIZ data are shown in the Supplement S2, S4). In general, there is good agreement amongst the methods and the magnitude of these differences is within the average within-method error (E w ), but the differences are important to highlight as they contribute to the between-method error (E b ). In the SAZ (Fig. 2a), the SOM-FFN differs from the other methods for summer and autumn from 1998 to 2008. The SVR method overestimates the seasonal amplitude pCO 2 (where the seasonal amplitude is the difference between the winter maxima and summer minima of pCO 2 ) relative to the other Table 2. A regional summary of the errors for the different methods. Note that the propagated errors are calculated as shown in Eq. (3) where the measurement and gridding errors are assumed to be constant at 5 µatm each (Pfeil et al., 2013;. The within-method and between-method errors are calculated using Eqs. (4) and (5)  methods for 2012 to 2014. In the PFZ (Fig. 2b), the SVR overestimates pCO 2 relative to the other methods during winter from 1998 to 2004, likely due to the method's sensitivity to sparse winter data . Figure 2c shows the time evolution of between-method errors for each biome. This panel highlights the seasonality of the estimates, specifically the increased heterogeneity of pCO 2 in summer and the impact this has on pCO 2 estimates. This is due to the more complex competing processes affecting pCO 2 during summer. To gain a better understanding of the seasonal processes, we consider the mean state of each season to characterise the drivers of opposing fluxes.

Ensemble seasonal cycle
The seasonal cycle of the pCO 2 for each biome (Figs. 2a-b and 3a-d) is coherent with seasonal processes reported in the literature (Metzl et al., 2006;Thomalla et al., 2011;Lenton et al., 2012. In all biomes, uptake of CO 2 is stronger during summer than in winter, giving rise to the strong seasonal cycle. This is due to the opposing influences of the dominant winter and summer drivers, partially damped by the seasonal cycle of temperature (Takahashi et al., 2002;Thomalla et al., 2011;Lenton et al., 2013). In winter, the dominant processes of mixing and entrainment results in increased surface pCO 2 and thus outgassing (Takahashi et al., 2009;Lenton et al., 2013;Rodgers et al., 2014). In summer, stratification allows for increased biological production and the consequent uptake of CO 2 , thus reducing the entrained winter DIC and associated pCO 2 (Bakker et al., 2008;Thomalla et al., 2011). However, stratification typically limits entrainment other than during periods of intense mixing driven by storms. This has an impact on primary productivity, DIC and pCO 2 (Lévy et al., 2012;Monteiro et al., 2015;Nicholson et al., 2016;Whitt et al., 2017).
The SAZ (Fig. 2a) is a continuous sink, where summer uptake (Fig. 3a) is enhanced by biological production and winter (Fig. 3c) mixing results in a weaker sink (Metzl et al., 2006;Lenton et al., 2012Lenton et al., , 2013. The same processes produce a similar seasonal amplitude in the PFZ (Fig. 2b), but stronger upwelling and weaker biological uptake result in a positive shift of the mean. This results in an opposing net summer sink and winter source. However, this is according to the mean state in the PFZ, and winter estimates of pCO 2 do in fact approach 0 µatm toward the end of the time series (Fig. 2b).
Also apparent from Fig. 3 is that, over and above the latitudinal gradient, pCO 2 is zonally asymmetric within each biome during summer (Fig. 3a), when biological uptake of CO 2 increases. Zonal integration of pCO 2 could thus dampen magnitudes of regional pCO 2 . A regional approach is therefore needed to examine the regional characteristics of seasonal and interannual variability of pCO 2 and to understand its drivers.
3.3 Regional pCO 2 variability: zonal and basin contrasts In Fig. 4, pCO 2 is decomposed into nine domains by biome and basin, with the boundaries defined in Fig. 1 (showing the SAZ and PFZ; air-sea CO 2 fluxes displayed in Fig. S4 in the Supplement). The regional estimates are plotted as time series for pCO 2 (black lines). The blue and orange lines show the respective annual maxima (typically winter) and minima (typically summer). The projected summer minima (dashed blue lines) are calculated by subtracting the mean seasonal amplitude from the winter maxima (Fig. 4). The projected summer minima are the expected summer pCO 2 under the assumption that summer pCO 2 is dependent on, but not restricted to, the baseline set by winter.
As found by Landschützer et al. (2015), the estimates of pCO 2 (Fig. 4) show that the Southern Ocean sink strengthens from 2002 to 2011 in all domains, referred to as the reinvigoration. This is preceded by a period of a net weakening sink (Fig. 4b, d, e) in the 1990s, referred to as the saturation period by Le Quéré et al. (2007). In other words, the ensemble shows the same trend found in past literature Ritter et al., 2017), but as with these studies we also find that there is large uncertainty in the interannual variability of the ensemble estimate, as shown by the between- c) the standard deviation between ensemble members for the three biomes, which is analogous to the between-method error (Eq. 5). The within-method (E w ) and between-method (E b ) errors are shown for each biome. For a more detailed breakdown of the errors see Table 2. method error in Fig. 4 (see Fig. S4 for the spread of the product estimates, including the Jena mixed layer scheme by Rödenbeck et al., 2014). This disagreement between methods is likely driven by the sparse coverage of pCO 2 measurements in the Southern Ocean, with empirical methods interpolating the sparse data differently Ritter et al., 2017). We thus present our methods and results as a framework to assess the drivers of interannual variability of pCO 2 . Key to understanding the mean interannual variability is that it is the net effect of the decoupled seasonal modes of variability for summer and winter. This is particularly evident in the PFZ (Fig. 4d-f). Here, and in the other biomes, the net strengthening of the CO 2 sink is mainly linked to a reduction of pCO 2 in winter for the majority of the time series. This corresponds with the findings of Landschützer et al. (2016), who linked the reinvigoration to the decadal variability of the SAM -the dominant mode of atmospheric variability in the Southern Hemisphere (Marshall, 2003). In contrast, summer pCO 2 variability is shorter (roughly 4-6 years), thus providing interannual modulation of longer-timescale winter variability. This is demonstrated well in the Indian sector of the PFZ, where a decrease in winter pCO 2 from 2002 to 2011 is offset by weakening of the summer sink from 2006 to 2010 (Fig. 4d). Similarly, in the Atlantic and Pacific sectors of the PFZ, decoupling occurs from ∼ 2011 to the end of 2014, with a rapid increase in the strength of the summer sink.
The mean amplitude of the seasonal cycle of pCO 2 -the mean difference between the summer minima and the win-  ter maxima -is perhaps a better way of understanding the strength of the seasonal drivers than the mean pCO 2 . For example, the Atlantic sectors of the SAZ and PFZ (Fig. 4c, f) have the strongest seasonal variability (14.11 and 25.83 µatm respectively). This contrasts with the relatively weak seasonal amplitude in the Indian sector of the Southern Ocean, which has mean amplitudes of 7.06 and 13.64 µatm for the SAZ and PFZ respectively (Fig. 4b, e). This contrast can also be seen by comparing the mean seasonal maps of pCO 2 in Fig. 3a and c. In summer, strong uptake in the eastern Atlantic sector of the Southern Ocean is indicative of large biological drawdown of CO 2 by phytoplankton (Thomalla et al., 2011). Conversely, relatively low primary production in the Indian sectors of the SAZ and PFZ result in a small seasonal amplitude (Thomalla et al., 2011). This large discrepancy in biological primary production is related to the availability of iron, a micronutrient required for photosynthesis. The lack of large land masses, which are a source of iron, in the Indian sector of the Southern Ocean could be a contributing factor to the lack of biomass (Boyd and Ellwood, 2010; Thomalla et al., 2011). Figure 4 gives us insight into the magnitude of interannual pCO 2 variability as well as the character of these changes, i.e. decoupling of interannual winter and summer modes of variability. This alludes to the point that pCO 2 is responding to different adjustments of seasonal large-scale atmospheric forcing and/or responses of internal ocean dynamics in the Southern Ocean (Landschützer et al., , 2016De-Vries et al., 2017).

Framework: seasonal deconstruction of interannual variability
In order to capture the decoupled 4-6-year short-term variability observed in summer, the estimates are divided into four objectively selected periods (P1 to P4). The periods are each 4 years long with the exception of P1, which is 5 years long due to the fact that the duration of the time series is not divisible by 4 (with a total duration of 17 years). Given a longer time series, this analysis would benefit from testing different durations for each period, as well as varying the starting and end years.
These four periods are too short for trend analyses (Fay et al., 2014), but the intention here is to identify periods that are short enough to resolve interannual changes of large-scale drivers of the winter and summer pCO 2 that would otherwise be averaged out over longer periods. We then calculate the relative anomaly between each successive period rather than an anomaly of the mean state (e.g. P2-P1). As a result, four periods give rise to three sub-decadal-scale transition anomalies for summer and winter: A (P2-P1), B (P3-P2), and C (P4-P3). We do this separately for each method rather than using the ensemble mean (see Sect. 3.1.4 for calculations). The mean of the method anomalies for each transition is then taken. These anomalies are considered significant if the absolute estimate of the anomaly is larger than the standard deviation between the methods for each period.
Note that although only summer and winter anomalies are discussed, it is recognised that autumn and spring could be equally mechanistically important. Winter anomalies of pCO 2 , wind stress, SST, and MLD are shown in Fig. 6, while summer anomalies of pCO 2 , wind stress, SST, and Chl a are shown in Fig. 7. Chl a is potentially a more important driver in summer than the generally shallow summer MLD (the omitted plots are shown in Figs. S5 and S6).

Uncertainty of transition anomalies
The transition anomalies are not calculated from the mean of the three products. Rather, we calculate the anomalies for each individual product with a n(p ) = S n(p) − S n(p−1) , where s are the estimates for a particular product, n represents an individual product and p represents P1 to P4. The result, a n(p ) thus represents the anomaly for two periods for a particular product. We then calculate the average of the anomalies with a p = 1 N · N n=1 a n(p ) , where N is 3, the number of products. We then calculate the standard deviation of the three anomalies (e p ), which is analogous to the between-method error, with e p = 1 N · N n=1 a n(p ) − a n(p ) 2 , where the terms are consistent with those above. We use e p as an uncertainty threshold where anomalies are only considered significant if |a p | > e p . These regions are masked in Figs. 6a-c and 7a-c. Figure 5 shows the winter (a-c) and summer (d-f) e p for each transition anomaly. Figure 6 shows the three transition anomalies for the four periods shown in Fig. 4. It is clear that there is large uncertainty around the pCO 2 anomalies in the Southern Ocean owing to the differences between the three empirical methods, caused by a paucity of in situ measurements of pCO 2 . However, there are still small regions that show anomalies with confidence. Figure 6d-l also show the Pearson's correlation coefficients for each of the driver variables with pCO 2 for the regions that are above the uncertainty threshold. The correlations in Fig. 6 show that MLD is a dominant predictor of pCO 2 in winter (Fig. 6j-l), with wind stress being a stronger predictor only in Transition B (6e). However, these correlations are all less than |0.3|, indicating that the relationship between pCO 2 and MLD is complex and nonlinear. Moreover, spatial inconsistency in the relationship between pCO 2 and the drivers reduce the correlations, which are applied for the entire domain (above the threshold). This is likely due to MLD being a metric that measures the complex interaction of heat, stratification, and mixing processes (Abernathy et al., 2011) -mechanisms relating to SST and wind stress. We now discuss the results by transition.

Drivers of winter pCO 2 variability
In Transition A (first column of Fig. 6), MLD is the strongest driver. Deeper mixed layers in the Pacific and eastern Indian sectors of the Southern Ocean correspond with increased deepening, correlating with increased pCO 2 . The reduction of pCO 2 along the boundary of the Atlantic and Indian sectors of the SAZ corresponds with increased SST. This agrees with the hypothesis put forward by Landschützer et al. (2015) that warmer SSTs in the Atlantic led to increased uptake of CO 2 . However, the same is not true for the western Indian sector of the PFZ, where cooling and deepening MLD results in a reduction of pCO 2 .
Increased uptake of pCO 2 across the boundary of the Atlantic and Indian sectors of the SAZ continues into Transition B (second column of Fig. 6). This is again accompanied by an increase in SST (Fig. 6h). The reduction of pCO 2 extends to the eastern Indian sector of the SAZ and Tasman Sea. This corresponds with weak shoaling of the MLD, weak warming and a reduction of wind stress (Fig. 6e, h, k). Conversely, in the eastern Pacific, cooling surface temperatures, weaker winds, and shallower MLDs correspond with a reduction of pCO 2 , again in agreement with . The large reduction of pCO 2 in the Indian sector of the PFZ corresponds with an increase in temperature; however, there is also an increase in the depth of the MLD -this interaction is mechanistically unlikely and may be an artefact of the sparse data in this region.
In Transition C, the reduction of pCO 2 in the Indian and western Pacific sector of the SAZ corresponds with warmer SST and shallower MLDs. Once again there is a region in the Indian sector of the PFZ that experiences a potentially spurious reduction of pCO 2 corresponding with deeper MLDs. The anomalies in the rest of the domain are not significant.

MLD-driven interannual variability of pCO 2 in winter
Our results indicate that there is not one dominant driver of pCO 2 interannual transition anomalies in winter. While MLD is on average the stronger driver, its dominance of pCO 2 is only marginal over SST and wind stress (Fig. 6). This marginal dominance over the two other drivers is likely due to MLD being a metric that integrates the complex interaction between wind-driven mixing and winter heat loss to the atmosphere, of which SST is a response (De Boyer Montégut et al., 2004;Sallée et al., 2010). Mechanistically, deeper MLDs would result in greater entrainment of DIC-rich deep waters, while shallower MLDs entrain less DIC-rich waters, thus reducing the DIC pool in winter resulting in potentially stronger pCO 2 uptake in the surface ocean .
An important point to note is that SST is negatively correlated with pCO 2 in Fig. 6g-i. This is contrary to what is expected for solubility-driven changes of pCO 2 (Takahashi et al., 1993). This indicates that SST -a response to underlying variability and trends in winter buoyancy and mixing -is not a driver of pCO 2 changes in most regions of the Southern Ocean. There are some small sub-regions, where SST could drive the pCO 2 trend, such as in the eastern Pacific sector of the PFZ during Transition B (Fig. 6b, h) but they are spatially and temporally limited. Our results suggest that, in winter, a complex interaction of changing wind stress and buoyancy fluxes that influence MLD and entrainment may play a stronger role than thermodynamics in explaining the pCO 2 interannual transitions.
Wind-stress anomalies (Fig. 6d-f) do not correlate strongly with pCO 2 anomalies, with the exception of Transition B, when it has the strongest correlation. We propose that this lack of coherence between the two variables may be a result of two compounding points. Firstly, wind stress is the only truly independent driver in the analysis, with SST and MLD both being used as proxies for pCO 2 in each of the products. Secondly, the wind stress shown in Fig. 6d-f considers only wind strength, so it does not take into account potential meridional changes in atmospheric circulation. This is the primary hypothesis presented in Landschützer et al. (2015), suggesting that atmospheric circulation became more zonally asymmetric. This induced a southward shift of warmer waters over the Atlantic and Indian sectors, reducing the depth of the MLD. Conversely, in the eastern Pacific cold winds induced colder SST and thus an increase in solubility.
Past studies have related the variability of Southern Ocean wind stress to the SAM, where the multi-decadal increasing trend has been cited as a reason for the saturation in the 1990s (Marshall, 2003;Le Quéré et al., 2007;Lenton and Matear, 2007;Lovenduski et al., 2008). The SAM is often represented as a zonally integrating index (Marshall, 2003), but more recent studies have shown that the SAM, as the first empirical mode of atmospheric variability, is zonally asymmetric (Fogt et al., 2012). The zonal asymmetry of the SAM is thought to be linked with the El Niño-Southern Oscillation and is strongest in winter, particularly over the Pacific sector of the Southern Ocean during a positive phase, in accord with the dipole nature of the Pacific-Indian winter windstress transition observed in Fig. 6d, e (Barnes and Hartmann, 2010;Fogt et al., 2012). Fogt et al. (2012) noted that the SAM has become more zonally symmetric in summer since the 1980s, matching the characteristics of the anomalies of wind-stress transitions seen in Fig. 6d-f and the hypothesis of Landschützer et al. (2015).
In summary, our analysis of the drivers of pCO 2 is consistent with the atmospheric asymmetry (dipole) conceptual model associated with the SAM proposed by Landschützer et al. (2015). However, our results suggest that interannual pCO 2 trends are explained by DIC dynamics rather than by the thermodynamic response of pCO 2 . A key part of this emphasis on DIC is that our results indicate that pCO 2 and SST are not correlated in a way that supports a thermodynamic control of pCO 2. The reasons for these differences are not clear at this stage, but they could include differences in the temporal resolution of the two studies: the resolution of the seasonal extremes in this study (seasonal modes) vs. annual mean in Landschützer et al. (2015).

Anomalies of pCO 2 and its summer drivers
Compared to the winter transitions, summer transitions (Fig. 7) have larger areas where the anomalies between products are within the bounds of the uncertainty. This may be due to the larger magnitude of the anomalies in summer compared to winter (Fig. 6). In summer we also see that Chl a (Fig. 7j-l) is likely the first-order driver with the highest correlation scores for transitions A and B.
Transition A (P2-P1 in Fig. 7a) is marked by a decrease of CO 2 in the SAZ (Tasman shelf region), coinciding with an increase in Chl a. The Drake Passage region experiences a strong reduction of pCO 2 in the PFZ, as found by Munro et al. (2015) and Landschützer et al. (2015). Unlike in the Tasman basin, this reduction of pCO 2 is not accompanied by a strong increase in Chl a, but rather a reduction of wind stress and an increase in SST. This is contrary to the annually integrated analysis of Landschützer et al. (2015), who found that cooling drove a reduction of pCO 2 in the eastern Pacific sector of the PFZ. This difference likely arises from the integration of seasons and a longer period (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) compared to the framework used in this study. In the Indian sector of the PFZ the products agree on a weak increase in pCO 2 corresponding with a weak reduction of Chl a. Transition B (P3-P2 in Fig. 7b) shows a large reduction of pCO 2 in the Atlantic sector of the PFZ and southern SAZ. The reduction coincides with an increase in Chl a and SST in the same region, in agreement with .
In Transition C (P4-P3 in Fig. 7c) the reduction of the pCO 2 is widespread in the Indian and Pacific oceans in both biomes, as the increase in Chl a is similarly widespread; however, the increase in Chl a in the Indian sector of the SAZ is not strong compared to other regions where pCO 2 and Chl a variability correspond. Conversely, there is a reduction in Chl a and concomitant increase in pCO 2 along the Polar Front in the Atlantic sector, coinciding with the position of the Antarctic Circumpolar Current (ACC) -a region with high eddy kinetic energy (EKE) (Meredith, 2016).
Based on these cases we suggest that pCO 2 is driven primarily by Chl a in regions with high-Chl a concentrations. Note that we will not try to explain Chl a variability, which is complex due to the multitude of factors influencing phytoplankton growth (Thomalla et al., 2011). We further suggest that in regions of low Chl a, buoyancy forcing and mixing are higher-order drivers. As suggested for winter variability, these two mechanisms are a complex interaction of variables of which SST and wind stress are a part (Abernathy et al., 2011).
This then raises the importance of the magnitudes of the interannual variability of pCO 2 and its drivers. For example, the anomalies of pCO 2 and SST are larger in summer than in winter. Conversely, the wind-stress anomalies are larger for winter than in summer. This is an important consideration for analyses that aim to understand the driving mechanisms, where annual averaging would weight seasonally asymmetric responses of pCO 2 and its drivers unequally. Figure 7. Relative anomalies of summer pCO 2 (a-c), wind stress (d-f), sea surface temperature (g-i), and mixed-layer depth (j-l) for four periods (as shown above each column). The thin black lines show the boundaries for each of the nine regions described by the biomes (Fay and McKinley, 2014a) and basin boundaries. Regions with dots are where pCO 2 anomalies are not significant; i.e. standard deviation of the anomalies between methods are greater than the absolute mean of method anomalies, as described in Eqs. (6)-(7).

Chlorophyll-dominated interannual anomalies of pCO 2 in summer
Our finding that Chl a is the dominant driver of interannual pCO 2 variability should not be surprising given that models and observations support this notion (Hoppema et al., 1999;Bakker et al., 2008;Mahadevan et al., 2011;Wang et al., 2012;Hauck et al., 2013Hauck et al., , 2015Shetye et al., 2015). However, our data show that the dominance of interannual Chl a variability over pCO 2 is largely limited to regions where Chl a is high, such as the Atlantic, the Agulhas retroflection and south of Australia and New Zealand (Fig. 8).  Thomalla et al. (2011) (using SeaWIFS data). Seasonality is calculated as the correlation between the mean annual seasonal cycle compared to the observed chlorophyll time series. A correlation threshold of 0.4 is applied to each time series to distinguish between regions of high and low seasonality; similarly, a threshold of 0.25 mg m −3 is used to distinguish between low or high chlorophyll waters. Black lines showing the fronts are calculated using altimetry thresholds from Swart et al. (2010).
The spatial variability of high Chl a regions in the Southern Ocean is complex due to the dynamics of light and iron limitation (Arrigo et al., 2008;Boyd and Ellwood, 2010;Thomalla et al., 2011;Tagliabue et al., 2014Tagliabue et al., , 2017. This complexity is highlighted in Thomalla et al. (2011), where the Chl a is characterized into regions of concentration and seasonal cycle reproducibility (SCR; Fig. 8). The SCR is calculated as the correlation between the mean annual seasonal cycle and the observed chlorophyll time series. Here we use the approach of Thomalla et al. (2011), in Fig. 8, as a conceptual framework to understand the interannual variability of pCO 2 .

High-chlorophyll regions
While regions of high SCR (dark green in Fig. 8) do not correspond with the interannual variability of Chl a (Fig. 7j-l), the framework by Thomalla et al. (2011) does present a hypothesis by which the variability of Chl a and its drivers can be interpreted, i.e. that the variability of Chl a in a region is a complex interaction of the response of the underlying physics (mixing vs. buoyancy forcing, which modulate light via MLD and iron supply) to the interannual variability in the drivers (SST and wind stress). This complexity is exemplified by strong warming in the Atlantic during Transition B, which results in both an increase and decrease in Chl a, with inverse consequences for pCO 2 . The effect is even stronger in Transition C, where strong cooling in the Atlantic results in both a decrease and increase of Chl a (Fig. 7i, l). In both transition A and B, the respective increase and decrease of Chl a occur roughly over the ACC, while the opposing effects during transitions A and B occur roughly to the north and south of the ACC region. These temperature changes may impact the stratification of the region, but complex interaction with the underlying physics results in variable changes in Chl a.
It is clear that, while there is a relationship between Chl a and pCO 2 as well as a relationship between wind stress and SST in summer, the relationship between wind forcing, Chl a, and pCO 2 is not as strong as in the winter anomalies (Fig. 6). It may be that enhanced summer buoyancy forcing resulting from summer warming and mixed layer eddies drives a more complex response to wind stress in the form of vertical velocities and mixing, which influence the iron supply and the depth of mixing (McGillicuddy, 2016;Mahadevan et al., 2012).
Mesoscale and sub-mesoscale processes may have a part to play in these dynamic responses of Chl a to changes in SST and wind stress (amongst other drivers). For example, eddy-driven slumping could act to shoal the mixed layer rapidly (Mahadevan et al., 2012;Swart et al., 2015b;du Plessis et al., 2017). This allows phytoplankton to remain within the euphotic zone, thus ensuring growth as long as iron is not limiting. Similarly, Nicholson et al. (2016) and Whitt et al. (2017) demonstrated that sub-mesoscale processes could supply iron to the mixed layer by sub-mesoscale mixing. Importantly, these mechanisms rely on a mixing transition layer that has sufficient iron to sustain growthweak dissolved iron gradients in the Pacific and eastern Indian sectors of the Southern Ocean could explain the lack of phytoplankton in these regions (Tagliabue et al., 2014;Nicholson et al., 2016). Much of the spatial character of the transition anomalies occurs at mesoscale, which strengthens the view that these mesoscale and sub-mesoscale processes may be key to explaining changes in Chl a (Fig. 7j-l).

Low-chlorophyll regions
Entrainment and stratification can explain much of the variability in the eastern Pacific and Indian sectors of the PFZ (with the exception of the wake of the Kerguelen Plateau). For example, in the eastern Pacific in Transition A (Fig. 7a,  d, g), strong warming and weaker winds have little impact on Chl a, but a decrease in pCO 2 is observed. Conversely, cooling in the western Indian sector of the PFZ results in a weak increase in pCO 2 during the same transition. In both these cases, the effect of cooling or warming on pCO 2 is negligible relative to the impact of entrainment or stratification respectively. The effect is reversed in the eastern Pacific during Transition B, where strong cooling results in a weak reduction of pCO 2 rather than the increase that would be expected from entrainment. This is the mechanism that Landschützer et al. (2015) ascribed to the reduction of pCO 2 in the Pacific, but the effect observed in Fig. 6b is weak.
In summary, regions with high-biomass Chl a integrate the complex interactions between SST, wind stress, MLD, and sub-mesoscale variability, resulting in large interannual pCO 2 variability compared to low-biomass regions, where wind-driven entrainment and stratification are more likely drivers of pCO 2 .

Synthesis
In this study, an ensemble mean of empirically estimated pCO 2 is used to investigate the trends and the drivers of these trends in the Southern Ocean. The estimated pCO 2 shows that the seasonal cycle is the dominant mode of variability imposed upon weaker interannual variability. The ensemble estimates are separated into domains defined by functional biomes and oceanic basins to account for the roughly basin-scale zonal asymmetry observed in preliminary analyses of pCO 2 (Fay and McKinley, 2014a). A seasonal framework is applied to the domains, revealing that winter and summer variability is decoupled for each region. The increase and subsequent decrease of pCO 2 (and air-sea CO 2 fluxes) is in accordance with recent studies showing a saturation of the Southern Ocean CO 2 sink in the 1990s, followed by the reinvigoration in the 2000s (Le Quéré et al., 2007;Landschützer et al., 2015).
While there is agreement around the mean of the ensemble, there is a large amount of uncertainty around the estimates due to a lack of agreement between products on a regional level. This uncertainty likely stems from the way that each method interpolates sparse winter data . We thus interpret only regions where the three empirical products are in agreement.
We suggest that changes in the characteristics of the seasonal cycle of the drivers of pCO 2 define the interannual variability of pCO 2 . In other words, the mechanisms that drive interannual modes of variability are embedded in the seasonal cycle.
Using this approach, we propose a refinement of the hypothesis put forward by Landschützer et al. (2015) by adding a seasonal constraint. The authors posit that pCO 2 variability is driven by changes in atmospheric circulation that in turn affect advection of water masses, thus impacting stratification. Our results also show that winter pCO 2 variability is best correlated with MLD, which indicates that en-trainment of deep DIC-rich water masses is an important mechanism of pCO 2 variability (Lenton et al., 2009Landschützer et al., 2015). The inverse relationship between SST and pCO 2 also suggests that in most cases pCO 2 is not thermodynamically controlled. Winter pCO 2 variability has a longer mode than summer variability, which we attribute to the decadal-mode variability of the Southern Annular Mode (Lovenduski et al., 2008;Fogt et al., 2012;Landschützer et al., 2016). This mechanism is likely dominant in winter due to its role in large seasonal net heat losses that drive convective overturning of the water column.
We suggest that interannual summer variability of pCO 2 occurs from a baseline set by an interannual winter trend. Moreover, the shorter-timescale summer interannual variability of pCO 2 (roughly 4-6 years) is driven primarily by Chl a. Buoyancy forcing and mixing still influence pCO 2 in summer but are lower-order drivers. We propose that the interannual variability of the summer seasonal peak is linked to the complex interaction of mid-latitude storms with the strong mesoscale and sub-mesoscale gradients in the Southern Ocean.
Overall, we propose that although mechanisms linked to winter wind stress explain the decadal trends in the strengthening and weakening of CO 2 uptake by the Southern Ocean, summer drivers may explain the shorter-term interannual variability (Lovenduski et al., 2008;Landschützer et al., 2015).
Lastly, it is important to note that this study can be improved by two factors. Firstly, increasing the length of the time series would allow for the identification of regular seasonal modes of variability. Moreover, the length of the anomaly periods could also then be adjusted to understand the variability of the drivers better. Secondly, improving machine-learning estimates of pCO 2 so that there is better regional agreement between products would decrease the area of insignificant variability. Rödenbeck et al. (2015) and Ritter et al. (2017) attribute the uncertainty to the individual methods' interpolation of sparse data in the Southern Ocean. This issue is being addressed by the community with autonomous sampling platforms closing the "observation gap". However, strategic deployment and sampling strategies will be critical to constrain and improve our understanding of CO 2 in the non-stationary context Monteiro et al., 2015).
Data availability. The data used in this study has been published  and can be accessed at https://doi.org/10.6084/ m9.figshare.5369038.v3 The Supplement related to this article is available online at https://doi.org/10.5194/bg-15-2361-2018-supplement.