Winter Air-Sea CO 2 Fluxes Constructed From Summer Observations of the Polar Southern Ocean Suggest Weak Outgassing

The Southern Ocean plays an important role in the global oceanic uptake of CO 2 . Estimates of the air-sea CO 2 flux are made using the partial pressure of CO 2 at the sea surface ( ), but winter observations of the region historically have been sparse, with almost no coverage in the Pacific or Indian ocean sectors south of the Polar front in the period 2004–2017. Here, we use summertime observations of relevant properties in this region to identify subsurface waters that were last in contact with the atmosphere in the preceding winter, and then reconstruct “pseudo observations” of the wintertime

. Observations of surf 2 pCO or its driving variables from both seasons are therefore required in order to constrain the air-sea flux.
Estimates of the air-sea flux of CO 2 from observations have traditionally used the sea surface partial pressure of CO 2 from ships and moorings. Existing observations of surf 2 pCO have been compiled in the Surface Ocean CO 2 Atlas (SOCAT, Bakker et al., 2016), and they are sparse in the Southern Ocean, particularly so in the wintertime and south of the Polar Front (see Figure 1). The available observations are used to generate time-evolving synoptic maps of surf 2 pCO by one of a number of gap-filling techniques (e.g., Gregor & Gruber, 2021;Landschützer, Gruber, Bakker, & Schuster, 2014;Rödenbeck, Keeling, et al., 2013;Takahashi, Sutherland, Wanninkhof, et al., 2009). The maps are then used in conjunction with atmospheric CO 2 data and a gas transfer parameterization to obtain the air-sea flux (see Section 4).
The sparsity of observations has meant that different observationally based estimates of the atmosphere-ocean CO 2 flux have their largest levels of disagreement in the Southern Ocean (Gregor et al., 2018;Rödenbeck, Bakker, et al., 2015). Despite this, an attempt by Gregor et al. (2017) to ascertain the effect of undersampling using model output found only a weak bias in ΔpCO 2 . However the applicability of this result to the real ocean depends on the representation of the Southern Ocean in models, which has been shown to be poor (Lenton et al., 2013;Mongwe, Chang, & Monteiro, 2016). Mongwe, Vichi, and Monteiro (2018) found that CMIP5 Earth system models show excessive ingassing south of 65°S compared to observations, and that models mostly disagree with one another and with the observations on the seasonal cycle of the CO 2 flux. It has therefore been suggested that more observations are needed to reduce uncertainties in the Southern Ocean carbon sink and its variability, in particular south of the Polar Front (Fay et al., 2014;Mongwe, Vichi, & Monteiro, 2018).
In an attempt to address the sparsity of observations in the Southern Ocean, autonomous profiling floats began to be deployed  (SOCCOM) project . The floats have pH sensors from which estimates of surf 2 pCO can be made, thereby substantially improving the spatiotemporal coverage of the Southern Ocean, especially in winter when the shipboard observations are most sparse. Using the SOCCOM float data, Gray et al. (2018) found stronger winter outgassing compared with previous estimates from shipboard surf 2 pCO , particularly between the Polar Front and the seasonal ice zone. Bushinsky, Landschützer, et al. (2019) (hereafter B2019) then combined the shipboard SO-CAT data with SOCCOM float data, and used the pCO 2 mapping methods of Landschützer, Gruber, Bakker, Schuster, et al. (2013) and Rödenbeck, Keeling, et al. (2013) (hereafter R2013) to estimate the flux. They found weaker outgassing than Gray et al. (2018), but still a significant reduction in the overall sink compared with estimates using only SOCAT data. They also obtained different results depending on whether they used a combination of SOCAT and SOCCOM data or SOCCOM data only, raising the question of whether a lack of agreement in the two data sets was due to different sampling or a possible bias. They tested the effect of a possible bias of 4 μatm in the float-derived pCO 2 by subtracting a uniform value before carrying out the pCO 2 mapping, and found that this did not eliminate the observed reduction in the estimated Southern Ocean CO 2 sink. However, it may be that there are as yet undiscovered biases in the float data, which require ground-truthing and calibration against the more established shipboard measurements to ensure their accuracy. Furthermore, even assuming it is accurate, the float-derived surf 2 pCO can do nothing to improve historical coverage.
In this study, we present an alternative method of addressing the paucity of wintertime observations of pCO 2 in the Southern Ocean, based on the method of Nomura et al. (2014). By utilizing historical subsurface summertime observations from the Global Ocean Data Analysis Project (GLODAP, Olsen et al., 2019) v2.2019, we extrapolate new "pseudo observations" of wintertime . These pseudo observations add wintertime coverage in all sectors of the Southern Ocean. We then combine them with observations of surf 2 pCO from SOCAT version 2019 and use a multiple linear regression (MLR) to produce a gap-filled time-evolving map of pCO 2 from which we construct a new estimate of the air-sea CO 2 flux. In the next section, we describe how we have constructed and validated our pseudo observations, in Section 3 we describe the mapping, in Section 4 we explain our estimation of the air-sea flux, and in Sections 5 and 6 we present and discuss our results.

Method
GLODAP contains profiles of various carbonate system parameters with depth, including DIC, Total Alkalinity (TA), as well as temperature (T) and salinity (S). Such observations were used by McNeil et al. (2007)  and hence estimate the air-sea flux of CO 2 . They have also be used by Nomura et al. (2014) to produce extrapolated estimates of surf 2 pCO in the wintertime of 2005 using summertime observations from the seasonal ice zone in the Indian Sector of the Southern Ocean, using the principle we now describe.
During winter, typically a deep mixed layer develops while the surface of the ocean cools and gases exchange with the atmosphere (Figure 2). During spring and summer the surface layer warms and restratifies, forming a pycnocline above the lower part of the winter mixed layer that isolates the winter water from the atmosphere. In the Southern Ocean south of the Polar Front, a temperature minimum layer (TML) forms because the winter surface water is cooler than both the summer surface layer above and the deeper waters below (Tomczak & Liefrink, 2005). Taking summertime profiles of temperature with depth from GLODAP, we can therefore use this TML to identify waters that were last in contact with the atmosphere during the previous winter. We make the assumption that the TML water is representative of surface conditions in Sep-MACKAY AND WATSON 10.1029/2020JC016600 3 of 17

Depth
Temperature/DIC DIC summer DIC winter T summer T winter TML CO 2 , O 2 etc. exchange with atm Phytoplankton tember, regarding this as the end of winter before the surface layer begins to warm and restratify. We search for a TML in the GLODAP profiles in all non-winter months, defined as October to May, and we generate pseudo observations for all of these months. The vast majority of TML profiles (84%) are from December to March, and we have validated pseudo observations from these months in Section 2.2. We also find that excluding pseudo observations constructed from profiles in April to May, which have the longest gap since the previous winter, has a negligible impact on our air-sea CO 2 flux estimate (these pseudo observations are included in the results presented).
We identify the TML using the following procedure. Taking profiles of temperature, salinity, and pressure from each year of GLODAP, we remove observations where the pressure is >500 dbar, so as to limit our search for a TML to those waters that might reasonably have been part of a winter mixed layer. The summer mixed layer depth (MLD) is then identified as the depth at which the potential density increases above a given threshold, "MLDtol", compared to the surface, and we then seek a temperature minimum between this depth and 500 m, which we define as our TML. We find that an MLDtol of 0.03 kg m −3 gives MLDs that agree well with a MLD climatology from Holte et al. (2017) in terms of their mean seasonal cycle and the spatial pattern of the long term mean for our study region. We vary MLDtol between 0.02 and 0.04 kg m −3 in MonteCarlo ensembles which we use to explore the uncertainties on our air-sea CO 2 fluxes (see Section 4 and Table C1), but our results are not sensitive to this choice.
To check that we are correctly identifying the surface winter water in the TML, we compare the TML T and S values with September surface layer (top 10 m) T and S from a gridded Argo product (Roemmich & Gilson, 2009). We find the values are well correlated, with Pearson coefficients of 0.83 and 0.7, and root-mean-squared errors (RMSE) of 0.96°C and 0.16 PSU for T and S, respectively (note that inaccuracies in the gridded Argo product in poorly sampled regions, particularly in sea ice, contribute to the RMSE).
Having identified the TML, we obtain the concentration of DIC within that TML from GLODAP. This is then adjusted for the effects of biological activity that has occurred in the winter water since it was last in contact with the atmosphere. We make this adjustment by assuming a Redfield Ratio RR C:O between the change in DIC and the change in oxygen between winter and summer due to biological activity: where AOU is the apparent oxygen utilization of the TML water, also from GLODAP, which is assumed to be the change in oxygen since the water was last ventilated in wintertime. To account for wintertime surface oxygen undersaturation in sea ice-covered regions we apply a correction to AOU of the form   corr b ice

AOU
(1 aC )AOU, where C ice is the sea ice concentration from the NOAA/NSIDC product (Meier et al., 2017;Peng et al., 2013), and the values of "a" and "b" are established during the validation with SO-CAT data described in the next section. We vary RR C:O between two possibilities, 106 138 and 106 150 , which are classic and revised Redfield Ratios quoted in Sarmiento and Gruber (2013), in our MonteCarlo ensembles. We have ignored any spatial or temporal variability in RR C:O .
In the previous paragraph, we have described a correction to AOU in the sea ice zone; however SOCCOM float data suggest that surface waters are undersaturated in oxygen over much of our study region in winter (Bushinsky, Gray, et al., 2017). Based on our analysis of the same float data, we apply a uniform correction of −13.5 μmol kg −1 to AOU used in the construction of pseudo observations in the region between the Polar Front and the sea ice. We find in the next section that this correction has a detrimental impact on our validation; however, we have applied it on the grounds that it is likely needed to remove a known bias. It may be that a spatially (and perhaps temporally) varying AOU correction would reconcile the discrepancies, but this would only be feasible for the most recent years since the deployment of the biogeochemical floats, and we will leave it for future work.
In order to create pseudo observations of trate (Ni) from the GLODAP gridded climatology. We use dissociation constants of Millero et al. (2006), Dickson (1990), and Lee et al. (2010), and finally combine the DIC, TA, Argo T and S, and GLODAP Si and phosphate (P) to calculate pCO 2 using CO2SYS. We take a mean over a surface layer, varied between 5 and 20 m in MonteCarlo ensembles, of the DIC, TA, T and S input to CO2SYS to produce our pseudo observations of surf 2 pCO ; for Si and P we have used fixed values (see Table C1). In converting our pseudo observations of DIC to surf 2 pCO in this manner, and in applying the AOU correction due to sea ice detailed above, we have assumed that our pseudo DIC from a summertime TML at a depth of between ∼50 and ∼200 m (see Figure A1) is representative of the surface DIC in the previous September within the same 1° latitude by 1° longitude box. This assumption excludes the possibility of horizontal transports and mixing having altered the properties of the TML over time, and we are not able to correct for this; however the validation we show in the next section gives us confidence that these effects are not important.

Validation
Where possible, we have validated our wintertime pseudo observations of DIC using real wintertime observations, also from GLODAP. For this validation, there is only one year (1998)  This estimate does not account for the possibility of spatial variations in the surface DIC trend, which may be weak or even negative locally. We compare annual means ( Figure 3c) rather than individual observations because the real and pseudo observations are not collocated; however, there is good agreement at a more granular level as well, as can be seen on Figure 3b. On the means, the difference between the pseudo and observed DIC is 7.1 μmol kg −1 for 1998 and 8.5 μmol kg −1 for 2005/2009. If we remove the −13.5 μmol kg −1 bias correction to AOU discussed in Section 2, the differences in the means are −3.3 μmol kg −1 for 1998 and 5.7 μmol kg −1 for 2005/2009. We additionally validate our pseudo observations with collocated, contemporaneous observations from SO-CAT. We identify SOCAT observations from the same month and year as a pseudo observation applies to, and which are located within ±1° latitude or longitude, and take a mean over the SOCAT observations to compare with that pseudo observation. We convert the SOCAT pCO 2 values to DIC and compare to our pseudo observations of DIC, thereby removing the effects of T and S on pCO 2 in the comparison. The conversion is done with CO2SYS, using inputs of T and S from SOCAT, assuming the same GLODAP climatological values of O, Si and Ni as were used to construct the pseudo pCO 2 , and estimating alkalinity using LIAR applied to those T, S, O, Si, and Ni (the climatological uncertainties on the latter three variables have negligible effect on this calculation). We convert each SOCAT observation collocated with a pseudo observation to DIC before taking a mean.
For the whole study period, there are 17 pseudo observations, from the years 2005 and 2008, for which there exist SOCAT observations for validation, and they are all in Drake Passage (Figure 4a). We find that overall agreement is good, with the pseudo observations slightly overestimating DIC compared with the 2005 SO-CAT observations (Figure 4b), and slightly underestimating it compared with those from 2008 ( Figure 4c). The RMSE and bias for both years combined are 13.9 and 4.1 μmol kg −1 , respectively. This reduces to 12.3 and 1.1 μmol kg −1 if we remove the AOU bias correction away from the sea ice. If we compare pCO 2 instead of DIC ( Figure B1), the RMSE and bias are 39.9 and 10.8 μatm, respectively. These reduce to 34.1 and 5.5 μatm without the AOU bias correction. While the bias correction has worsened the agreement between our pseudo observations and both the GLODAP and SOCAT observations we have used for validation, we will account for this by covering the range of AOU uncertainty in our MonteCarlo ensembles (see Section 4).   ed compared with SOCAT. These stations are in the region between the sea ice and the Polar Front where a bias correction has been applied to the AOU involved in the pseudo observations calculation; with this correction removed the four most northerly SOCAT observations fall within the error bars of the pseudo observations. Two further pseudo observations in 2005, one at around 60.7°S and one at 62.8°S, stand out as overestimates compared with SOCAT (the discrepancy is more obvious in the pCO 2 comparison shown on Figure B1). These pseudo observations derive from summer profiles where the temperature minimum was at 150 and 200 dbar, much deeper than the ∼70 dbar for the majority of the TML profiles in that year (2005, see Figure A1). However, the temperature profiles do still have a characteristic TML shape. We suggest that these deep temperature minima result from localized, short-lived deep winter mixing events which entrain deeper waters rich in natural carbon and thus produce pseudo observations with high values of DIC and pCO 2 . If these events did not occur at the right time to coincide with the validating SOCAT observations, this would explain the overestimate.
We use the comparison with SOCAT to establish optimized values of the parameters "a" and "b" used to correct the AOU in sea ice-covered regions. We explored a range of values for both parameters, creating sets of pseudo observations and calculating the RMSE for all 17 stations on Figure 4 for each set. The parameters that give the best fit are a = 0.9 and b = 0.3, with RMSE = 13.9 μmol kg −1 . The same optimization carried out using the pCO 2 values shown in Figure B1 results in the same best fit parameters, and an RMSE = 39.9 μatm.

pCO 2 Mapping
Having created pseudo observations that improve wintertime coverage of observed pCO 2 in the Southern Ocean south of the Polar Front (Figure 1), the next step is to use the available observations to map surf 2 pCO in space and time, as with the studies referred to in Section 1. We use a simple MLR, which while unable to fit the observations as closely as more complex methods such as the neural network approach of Landschützer, Gruber, Bakker, Schuster, et al. (2013), has the advantages of simplicity and modest computation requirements, allowing us to carry out our MonteCarlo ensembles to explore uncertainty. The variables on which we base the MLR are the atmospheric CO 2 mole fraction atm 2 xCO , T, S, and MLD. depends directly on T and S, and indirectly on MLD because surface mixing entrains DIC into the mixed layer from below. Including atm 2 xCO allows us to capture some of the interannual variability and the trend, however, as we will see later the former tends to be underestimated by this method. T and S have the advantage, thanks to the Argo program, that they are comparatively well known in time and space.
We use the Roemmich-Gilson (R-G) Argo product in which T and S are gridded with 1° horizontal resolution, on 58 depth levels (34 in the top 500 m), and with monthly values for our study period. We obtain atm 2 xCO from the NOAA marine boundary layer product (Dlugokencky et al., 2017), also with monthly values on a 1° × 1° grid. We obtain MLDs from the monthly climatology of Holte et al. (2017), which applies an algorithm to density profiles based on Argo data, also on a 1° × 1° grid. The MLD algorithm has proven more accurate than threshold methods, which tend to overestimate MLDs particularly in polar regions (Holte & Talley, 2009). We apply the same climatology to every year, thereby neglecting any interannual variability in MLD. There are gaps in the climatology which we fill by trilinear interpolation of the fields (dimensions of latitude, longitude, and month). We also tested filling the gaps with MLD values derived from a density threshold method applied to the R-G Argo product gridded T and S fields, and this gives nearly identical results.
In addition to our pseudo observations, we calculate surf 2 pCO directly from GLODAP observations, and we also take it from a gridded version of SOCAT. The GLODAP surf 2 pCO is constructed by calculating profiles of pCO 2 using CO2SYS on GLODAP profiles of DIC, TA, T, and S, and using the same Si and P as for the pseudo observations, and averaging those pCO 2 profiles over the 10 m surface layer. All three sets of observations (pseudo, GLODAP, and SOCAT) are combined and regressed on their associated T, S, MLD and year, to produce a time-evolving map of derived directly from GLODAP observations, the T and S is taken from the GLODAP profiles, again averaged over the surface layer, and the MLD is obtained from the climatology. In the case of the SOCAT gridded surf 2 pCO , T and S comes from the SOCAT product and MLD from the climatology. In all cases atm 2 xCO is taken from the nearest latitude/longitude grid cell to the pseudo, GLODAP or SOCAT observation in the relevant month and year.
Finally, before calculating the MLR we split the observations into two regions, one north and one south of the mean position of the Polar Front. This is a simplification of the approach used by Landschützer, Gruber, Bakker, Schuster, et al. (2013) and others, where the ocean is split into biogeochemical provinces, or "biomes"; in our case since we have relatively few pseudo observations (760), if we were to use similar biomes there would be several that contained no pseudo observations. Therefore, we split only at the Polar Front on the grounds that the pseudo observations can only be calculated to its south. We carry out the MLR both with and without the pseudo observations included in the data set, so as to determine their impact. We validated the MLR mapping by comparing the mapped gridded with the nearest in situ observations from the entire SOCAT database in the same month and year, obtaining an RMSE of 24.5 μatm. These observations are not independent of the mapping because they are the basis on which the SOCAT gridded values we include in our MLR are constructed, but this validation gives an indication of how well we are fitting the data.

CO 2 Flux Calculation
We calculate the air-sea flux of CO 2 from the difference between In Equation 2, k is the gas transfer velocity calculated according to Nightingale et al. (2000) and using CCMP 6-hourly winds (Atlas et al., 2011) and sea surface temperature from the surface layer of the R-G gridded Argo product, and K 0 is the solubility calculated according to Weiss (1974) and using sea surface temperature and sea surface salinity from the gridded Argo product.
atm 2 pCO was calculated from atm 2 xCO following Cooper et al. (1998) and using gridded Argo sea surface temperature and NCEP/NCAR sea level pressure (Kalnay et al., 1996).
Uncertainty on CO2 F comes from a number of sources, which we describe here, and summarize in Table C1. Each of the parameters and errors identified are varied in a 200-member MonteCarlo ensemble to estimate the envelope of uncertainty on CO2 F , and we also give the range of parameters explored in Table C1.
The GLODAP AOU used to adjust our pseudo DIC for biological activity between winter and summer may be overestimated because of wintertime undersaturation of surface oxygen (ΔO 2 ) in some regions. We assign normally distributed random errors to AOU with a mean of −13.5 μmol kg −1 and a standard deviation of 15 μmol kg −1 away from sea ice, and normally distributed random errors with a mean of zero and a standard deviation of 15 μmol kg −1 in sea ice (where the AOU correction described in Section 2 based on sea ice concentration has also been applied). The normally distributed errors are designed to cover the range of likely ΔO 2 values that would necessitate a correction to AOU. Between the Polar Front and the sea ice zone the ±1σ range is between −28.5 and +1.5 μmol kg −1 , and in the sea ice zone it is between ∼−56 and ∼−26 μmol kg −1 (the latter range is a crude estimate based on the effect of our sea ice-based correction on typical TML AOU values of ∼54 μmol kg −1 found in the GLODAP data and mean September sea ice concentration of 0.56 within the region of >0 sea ice concentration).
Values of T and S from the R-G gridded Argo product are used in multiple parts of the method: in the construction of pseudo observations through (1)  ; and in the calculation of (4) k and (5) K 0 used in estimating the air-sea CO 2 flux. Throughout we generate normally distributed random errors on T and S with MACKAY AND WATSON 10.1029/2020JC016600 8 of 17 a mean of zero and allow them to propagate, but in (1 and 2) we use standard deviations for the normal distributions based on the comparison of TML with surface Argo T and S from Section 2, whereas in (3-5) we use more modest estimates of the uncertainty (see Table C1).
The MLDs from Holte et al. (2017) have relatively large uncertainties and contribute significantly to the uncertainty on the flux. At each grid point in the MLD climatology, we use the standard error of the MLDs contributing to the climatological mean at that point as the standard deviation for a normal distribution of random errors applied to the MLD.
We vary between 5 and 20 m the depth of the surface layer over which we average the various quantities used to calculate surf 2 pCO in our MonteCarlo ensemble. We establish bounds on the values of silicate and phosphate input to CO2SYS by taking a mean over our study region of the upper and lower bounds of those parameters in the GLODAP climatology, and vary them between their upper, mean, and lower bounds. We also investigated the impact of varying the GLODAP climatological silicate, oxygen, and nitrate fields input to LIAR for estimating TA for the pseudo observations. Varying these inputs within the climatological uncertainties has a negligible effect on our results, and they are kept constant in our ensembles. Finally, we generate a normal distribution of multipliers for the gas transfer coefficient k to simulate an uncertainty of 20% on this parameter.

Results
South of the Polar Front and north of the seasonal ice zone, the annual mean CO2 F is generally out of the ocean on the eastern side of Antarctica and into the ocean on the western side (Figures 5a and 5b). Within the seasonal ice zone the annual flux is generally into the ocean, but small. The addition of the pseudo observations has increased the outgassing in all sectors of the ACC, most strongly near the edge of the seasonal ice zone around 60°E (Figure 5c). In winter, there is outgassing around most of the ACC between the Polar Front and the seasonal ice zone, and very little flux within the seasonal ice zone as expected (Figures 5d  and 5e). The impact on the winter mean of including the pseudo observations follows a similar spatial pattern as for the annual mean, but is larger in magnitude (Figure 5f; note that the color scale in subplots c and f differs from the rest). Nomura et al. (2014) reported sea-air flux densities of −0.4 mol C yr −1 m −2 in the seasonal ice zone between 32°E-58°E and 64°S-67°S for the winter of 2005. At those longitudes our estimate shows negligible flux at the extreme southern edge of our domain, which just overlaps with the northern extent of theirs.
The 2004-2017 mean seasonal cycle of CO2 F integrated over the region south of the Polar Front with the inclusion of the pseudo observations shows stronger outgassing/weaker uptake compared with that without pseudo observations (Figure 6a). We find a peak outgassing in July of 0.14 Pg C yr −1 (lower bound, upper bound 0.11, 0.17 Pg C yr −1 ) with the pseudo observations, and 0.04 (0.03, 0.05) Pg C yr −1 without. The overall annual mean flux is −0.02 (−0.04, −0.00) Pg C yr −1 with the pseudo observations, and −0.09 (−0.10, −0.09) Pg C yr −1 without. Compared to flux estimates based on the surf 2 pCO products of R2013 (green dashed line), B2019 (dark red dashed line), and Landschützer, Gruber, and Bakker (2017) (bright red dashed line, hereafter L2017), our results exhibit a slightly weaker seasonal cycle, rather similar in shape to the L2017 estimate. Note that the flux estimates presented on Figure 6 from our methodology and from those of R2013, B2019, and L2017 are exactly comparable, in the sense that we have calculated CO2 F using the mapped  which used the interpolation method of Rödenbeck, Keeling, et al. (2013), also applied to SOCAT data. The dark red line is based on  Compared to the estimates based on R2013, B2019, and L2017, our estimate without pseudo observations follows the general trend between 2004 and ∼2013 but with less variability, while the estimate with pseudo observations has a steeper downward trend. After 2013 the estimates based on other studies trend upwards while our results continue to follow the downward trend. The B2019 estimate (dark red dashed line) includes the use of SOCCOM float data from 2014, and this also has an influence on previous years which diminishes going back in time. This accounts for some of the strong upward trend, but the R2013 and L2017 estimates are based solely on SOCAT data, and they also trend upwards, by contrast to our results. Part of this discrepancy can be explained by the differences in gap-filling methods: we have used an MLR, which is constrained toward linearity and cannot reproduce the sharp change in trend captured by the more sophisticated techniques. If we split our MLR into two halves, training it first on the data from 2004 to 2011 and then from 2011 to 2017, we see an upward trend with the pseudo observations from 2011, that is still weaker than the estimates based on L2017 and B2019, but that overlaps with R2013 in 2017 ( Figure D1). The upward trend in the later period is slightly stronger with the pseudo observations than without. In the earlier period, the outgassing is stronger, the difference due to the pseudo observations is more pronounced, and the downward trend is steeper compared with the estimates based on B2019 and L2017. In this period, our estimate without pseudo observations quite closely follows the trends of L2017 and B2019. We discuss the differences between our results and those based on comparable earlier studies further in the next section.

Discussion
We have found that the addition of wintertime pseudo observations of The downward trend (increased uptake) in the annual mean flux from 2004 to 2017 seen in our estimates on Figure 6b is expected given the steady increase in atmospheric CO 2 concentrations over the same period; however it is at odds with other estimates based on surf 2 pCO reconstructed using more sophisticated techniques, which suggest a reversal of the downward trend in the latter part of the period. Because our estimates are constrained toward linearity by the linear nature of the MLR, we have tested the impact of splitting the MLR in two at 2011. We find that this does reproduce some of the upward trend, but it is weaker at least than that of fluxes based on surf 2 pCO reconstructions of L2017 and B2019. These estimates both used the neural network technique of Landschützer, Gruber, Bakker, Schuster, et al. (2013) for the reconstructions, but the former is based soley on SOCAT surf 2 pCO data and the latter includes data from SOCCOM floats. Another potential limitation to our method is that we have used a MLD climatology, which means that our CO2 F estimate is not influenced by interannual variability or trends in the MLDs. However, we find no obvious trend in MLDs estimated from the R-G gridded Argo product that would explain the upward trend in CO2 F suggested by other products.
Although the wintertime spatial coverage is much improved by the pseudo observations (Figure 1), the total number of grid points (lat/lon/month/year) occupied is still small: 760 compared with 27,983 from SOCAT. We tested the impact on our CO2 F estimates of artificially inflating the number of pseudo observations by simply including multiple copies of each pseudo observation in the data used in the MLR. We calculated CO2 F with up to 30 copies of the pseudo observations, such that they would be as numerous as the SOCAT observations. We find that increasing the influence of the pseudo observations increases the strength of the downward trend in the flux, therefore we conclude that they have not been swamped by the more numerous direct observations of The largest difference was found between our estimates and one based on a neural network reconstruction of surf 2 pCO from B2019 including SOCCOM data, in the years 2015-2017 when the float data have the most influence. We attempt a comparison between our pseudo observations and SOCCOM float-derived DIC, and are able to locate three examples of pseudo observations in the same month and year and within ±2° latitude or longitude of a SOCCOM float with derived DIC observations. For these three, the DIC of the pseudo observations is on average only 2.0 μmol kg −1 lower than the comparable SOCCOM float values; however, this is a very small number of observations from which to draw inferences from the comparison. B2019 also investigated the effect of a possible bias in the float-derived surf 2 pCO by subtracting a uniform value of 4 μatm from their estimate, concluding that a mean bias of that magnitude would not account for the difference made by the float data to the CO 2 flux. However, it may be that such a uniform correction does not account for a bias that is unevenly distributed, or for a larger bias as yet unidentified. Our results reinforce the suggestion made by B2019 that more work may be needed to determine the accuracy of float-derived surf 2 pCO .
The seasonal cycle of our estimated CO 2 flux south of the Polar Front (Figure 6a), with its low in January and high in July, is very similar to the seasonal cycle of the mean MLD for the same region, suggesting the MLD is a key driver. By contrast, the seasonal cycle of the mean surface temperature is out of phase, having its low in September and high in February. The seasonal cycle of the surface salinity, which is linked to surf 2 pCO through its relationship to alkalinity, is almost in phase with the CO 2 flux, with a low in February and a high in September. Gregor et al. (2018) conclude that winter variability of surf 2 pCO in the Southern Ocean is driven by stratification and mixing, whereas summer variability is driven by primary production. They also find MLD to be the dominant predictor of winter , consistent with our results. This dependence on MLD has been linked to the entrainment of DIC-rich deep waters in wintertime (Lenton et al., 2013). During our validation we found that pseudo observations with anomalously high DIC were derived from summer profiles with particularly deep TMLs, suggesting that isolated deep winter mixing events might be important for observed winter outgassing.
We have presented a novel estimate of the air-sea flux of CO 2 in the Southern Ocean south of the Polar Front that utilizes summertime subsurface observations of carbonate system parameters to boost the wintertime coverage of surf 2 pCO data which are required to estimate the flux. We find that these additional pseudo observations of surf 2 pCO result in an increase in the estimated winter outgassing in the region, but that this increase is largest for the mid 2000s, and reduces with time. In 2017, the most recent year of our analysis, we estimate a CO 2 flux of −0.14 Pg C yr −1 (lower bound, upper bound −0.16, −0.12 Pg C yr −1 ), which increases to −0.08 (−0.10, −0.05) Pg C yr −1 when the MLR is carried out separately on 2011-2017 data. This compares with +0.06 Pg C yr −1 derived from surf 2 pCO constructed using a neural network technique applied to a combination of data from ships and autonomous profiling floats. The difference is comparatively small, but it is significant that improving winter coverage does not bring our estimate in line with one which used the float data, so the origin of the discrepancy between shipboard and float data in the Southern Ocean winter remains unsolved. This is important because placing constraints on the Southern Ocean air-sea CO 2 flux and determining its variability is critical to our understanding of the global carbon cycle, and hence to our ability to make accurate climate projections.