Storms drive outgassing of CO2 in the subpolar Southern Ocean

The subpolar Southern Ocean is a critical region where CO2 outgassing influences the global mean air-sea CO2 flux (FCO2). However, the processes controlling the outgassing remain elusive. We show, using a multi-glider dataset combining FCO2 and ocean turbulence, that the air-sea gradient of CO2 (∆pCO2) is modulated by synoptic storm-driven ocean variability (20 µatm, 1–10 days) through two processes. Ekman transport explains 60% of the variability, and entrainment drives strong episodic CO2 outgassing events of 2–4 mol m−2 yr−1. Extrapolation across the subpolar Southern Ocean using a process model shows how ocean fronts spatially modulate synoptic variability in ∆pCO2 (6 µatm2 average) and how spatial variations in stratification influence synoptic entrainment of deeper carbon into the mixed layer (3.5 mol m−2 yr−1 average). These results not only constrain aliased-driven uncertainties in FCO2 but also the effects of synoptic variability on slower seasonal or longer ocean physics-carbon dynamics. Storms dominate the subpolar Southern Ocean, where upwelling CO2 drives outgassing that impacts global CO2 budget, yet how storms modify this outgassing is unknown. Here, the authors present coupled atmosphere-ocean observations to show how storm-driven ocean mixing and circulation cause substantial CO2 variability and outgassing.

T he Southern Ocean is a key component of the Earth's carbon budget. It accounts for 40-50% of the total mean annual ocean uptake of anthropogenic CO 2 (~1 Pg C yr −1 ) [1][2][3][4] . In addition, an increase of the annual mean outgassing of natural CO 2 from the Southern Ocean in the 1990s and the subsequent decrease in outgassing in the 2000s (resulting in a~0.5 Pg C yr −1 reinvigoration of ocean uptake) showed that the global ocean carbon budget is sensitive to variability in the Southern Ocean [5][6][7][8][9][10] . This outgassing variability has been linked to climate-mode forced variations in wind-driven upwelling that comprises the surfacing of deep waters with high concentrations in dissolved inorganic carbon (DIC) resulting in widespread but variable outgassing of CO 2 in the subpolar region, which counteracts the CO 2 uptake flux 2,9,10 . However, it is not well understood what role the daily-to-seasonal physics of the surface mixed layer, the critical boundary between the atmosphere and the upwelled reservoir of DIC, plays in modulating the magnitude of this outgassing flux 11,12 .
The subpolar outgassing region, south of the Polar Front, coincides with the core of Southern Hemisphere storm tracks (~50-65°S) 13,14 and is regarded as the windiest region in the world. Strong storms occur in regular succession (4-8 days 15 ) throughout the year, inflicting intense (>0.8 N m −2 16 ) but short-lived (2-3 days 15 ) surface wind stress. In the mid-latitudes of both hemispheres, storms energize the mixed layer at synoptic timescales (1-10 days), triggering enhanced vertical mixing (entrainment) and advection that result in surface and subsurface ocean exchanges of heat, momentum and chemical properties. Thus, storms drive significant variability in biogeochemistry 11,[17][18][19][20] and the amplitude and timing of seasonal mixed-layer depths 21,22 . Yet, despite these well-documented upperocean responses, there remains limited understanding of what impact these frequently passing storms have individually or cumulatively on the outgassing magnitude and variability in this globally critical region, the subpolar Southern Ocean.
The lack of understanding of how storms impact this region, may be further reflected by the strong biases in the seasonal cycle of ocean CO 2 in both Earth System Models and forced ocean models in the Southern Ocean, which exclude sub-grid scale mixed-layer dynamics [23][24][25] . Model experiments suggest that the exclusion of processes associated with storms (e.g. near-inertial waves) significantly reduces vertical mixing and the surface ocean supply of DIC, reducing winter outgassing and increasing summer uptake of CO 2 11 . Moreover, recent high-resolution (10-day) observations from biogeochemical floats have revealed an underestimation of the magnitude of CO 2 outgassing in shipbased contemporary estimates due to significant spatial-temporal aliasing by these observations 26,27 . Autonomous surface vehicle sampling at higher resolution (hourly or faster) shows further temporal aliasing of floats (i.e. a 10-day sampling interval may result in~20% uncertainty in the mean CO 2 flux). This aliasing is likely a result of unresolved mixed-layer responses to highly variable synoptic events 17,28 . There remains a significant gap in the understanding of the mechanisms that drive variability on synoptic timescales and how this synoptic variability rectifies on the seasonal cycle and mean of CO 2 fluxes.
Here, we address these gaps through high-resolution atmosphere-ocean observations from autonomous vehicles in the Atlantic sector of the subpolar Southern Ocean. We find that storms modulate the direction and magnitude of the air-sea CO 2 gradient (ΔpCO 2 ) and flux (F CO2 ) through two mixed-layer processes. Synoptic variability in the wind-driven Ekman flow advects upwelled DIC-rich waters from the south (and low DIC waters from the north) and explains most of the ΔpCO 2 variability on timescales from a few days to months. In addition, intense wind events drive turbulent entrainment of the deeper DIC reservoir into the surface mixed layer and strong outgassing F CO2 . We construct a process model, which captures the observed ΔpCO 2 variability to estimate its relevance across the subpolar Southern Ocean and on longer timescales. We show that Ekmandriven synoptic variance of ΔpCO 2 is spatially modulated by large-scale ocean fronts, where meridional DIC gradients are strong hence Ekman flows drive large DIC transports. However, Ekman-driven synoptic variance of ΔpCO 2 is largely oscillatory with little additive effect on the mixed-layer DIC budget on timescales longer than a couple of months. On the other hand, storm-driven entrainment, which is spatially modulated by variations in seasonal stratification, adds up to be relevant for the annual mean mixed-layer DIC budget. These results support the hypothesis that storm-driven ocean physics is a significant participant in the Southern Ocean carbon cycle.

Results
High-resolution observations of mixing and CO 2 variability. The subpolar Southern Ocean is where the upwelling of Upper Circumpolar Deep Waters (UCDW) supplies elevated concentrations of DIC to the base of the mixed layer 29 (Fig. 1). This gives rise to a circumpolar zonal band of spatially and seasonally varying outgassing of natural CO 2 (Fig. 1a, b). A surface autonomous vehicle (Wave Glider measuring atmospheric weather and carbon chemistry) and an ocean profiling vehicle (Slocum glider measuring physical characteristics of the water column) sampled one site (54°S, 0°E) in a coordinated way to yield an integrated high-resolution view of the interactions between the air-sea fluxes and upper-ocean dynamics in the South Atlantic Ocean (Fig. 1, refer to Methods).
The air-sea flux of CO 2 is defined (Eq. 1) as the product of the air-sea gradient of pCO 2 (ΔpCO 2 = pCO 2sea − pCO 2atm , where pCO 2sea is the partial pressure of CO 2 in the surface ocean and pCO 2atm in the atmosphere), the temperature-dependent solubility K s and the wind dependent air-sea gas transfer velocity K w 30 .
Notably, K w which is quadratic in wind speed 30 and ΔpCO 2 are both crucial to F CO2 , but this study focuses primarily on how storms influence ΔpCO 2 . The direction and part of the magnitude of F CO2 are set by ΔpCO 2 , the variability of which is dominated by pCO 2sea in our observations ( Supplementary Fig. 1). The observed pCO 2sea and hence ΔpCO 2 varied by~±10 µatm, as F CO2 oscillated between uptake and outgassing on synoptic timescales (1-10 days) (Fig. 1e). Several of the outgassing events coincided with the passage of storms (Fig. 1e-compare grey and red shaded areas). To put these results into perspective, the synoptic variability of ΔpCO 2 (about 20 µatm from peak to trough) is similar in magnitude to the seasonal amplitude of ΔpCO 2 and pCO 2sea for the subpolar Southern Ocean 10,31 .
To elucidate the causes of synoptic variability of observed ΔpCO 2 and pCO 2sea , we begin by decomposing the drivers of changes in pCO 2sea into the relative contributions made by thermal (pCO 2-SST ) and by non-thermal (pCO 2-DIC ) components 32 referenced from the start of the deployment (Fig. 1f). This decomposition shows two timescales of variability, which are key to explaining the mechanisms: first, the synoptic-scale variability of ΔpCO 2 and pCO 2sea (Fig. 1e, f) is dominated by the pCO 2-DIC component, which we will show is primarily caused by winddriven DIC transport in the ocean. In fact, the synoptic pCO 2-DIC (pCO 2-DIC ′, computed by removing the 10-day rolling mean from pCO 2-DIC ) explains about 70% of the variations in pCO 2sea (r 2 = 0.71, and refer to Figs. 1e, 2d, and Supplementary Fig. 1). Second, both pCO 2-DIC and pCO 2-SST also show summer seasonal trends, which are weakening for pCO 2-DIC and strengthening for pCO 2-SST (Fig. 1f). Over the duration of the deployment, these trends in pCO 2-DIC and pCO 2-SST ultimately yield changes of ~20 µatm, which are comparable in magnitude to the synoptic variability (Fig. 1e). The increasing trend in pCO 2-SST is due to progressively increasing sea surface temperature (SST) linked to seasonally driven solar warming ( Supplementary Fig. 3a). Meanwhile, the weakening trend in pCO 2-DIC is likely a consequence of a gradual decrease in DIC by~10 mmol C m −3 (or equivalently 1000 mmol C m −2 assuming the top 100 m is mixed) over 2 months due to biological productivity. Consistent with this inferred DIC drawdown of 1000 mmol C m −2 , bio-optical estimates of net primary productivity from sensors on the Slocum glider and satellites range from~13-40 mmol C m −2 d −1 or equivalently 700-2300 mmol C m −2 over the 56-day deployment ( Supplementary Fig. 4). Ocean advection and a net freshwater flux due to precipitation could also influence this trend in pCO 2-DIC significantly ( Supplementary Fig. 3b), but the relative contributions of advection and freshwater fluxes are not investigated in this study (refer to 31 ). Regardless of the driving mechanisms, the seasonal trends in pCO 2-DIC and pCO 2-SST approximately compensate for each other, thus ΔpCO 2 fluctuates about its initial value near zero, and the sign of ΔpCO 2 changes on synoptic timescales during the 2-month glider deployment at this site.
To link the synoptic variability of ΔpCO 2 and pCO 2sea (Fig. 1e, f) with wind-driven mixed-layer processes that connect the surface ocean to the DIC-rich subsurface reservoir (Fig. 1c)  surface friction velocity due to the wind (u * defined in Methods) and the dissipation rate of turbulent kinetic energy (ε) observed by the profiling glider (Fig. 2a, b). The dissipation rate is a measure of turbulence and is related to turbulent diffusion and vertical mixing through stratification. The vertical extent of the elevated ε is representative of the turbulent boundary layer depth, hereafter referred to as the mixing-layer depth (XLD, defined in ref. 33 ). Turbulent, deep-reaching mixing events down to 150 m correspond with strong wind (high u * ) compared with periods of shallower (about 50 m), weaker mixing events occurring during weak winds (low u * ). Here, the observed variability in the magnitude (ε) and depth of mixing (XLD) was driven primarily by wind (u * ), r 2 = 0.76 and r 2 = 0.75, respectively (Supplementary Figs. 5 and 6). In contrast to the highly variable XLD, the density-derived mixed-layer depth (MLD, defined in Methods) was nearly constant through mid-January at around 130 m (Fig. 2b). This initial MLD of 130 m is set by the winter MLD maximum, and it is coupled to the upper bound of the high-salinity UCDW (σ= 27.3 kg m −3 , defined in ref. 34 ). During mid-January the MLD shoaled to about 75 m, decoupling from UCDW and its carbon-rich waters (Fig. 2b, c). The MLD was not directly sensitive to the variability of the wind (r 2 = 0, Fig. 2a, b) and was thus distinctly different from the XLD. This is consistent with our conventional density threshold definition of the MLD (refer to Methods), which is not expected to vary with the XLD on synoptic or shorter timescales 33 . The absence of correlation between the MLD and the wind is also consistent with the model of Whitt et al.
(2019) 22 , in which sub-seasonal u * is weakly correlated with subseasonal MLD both in the subpolar Atlantic sector of the Southern Ocean and globally. In particular, since the MLD is set by the integrated effect of mixing, it does not vary on the same timescales as the wind. Thus, the MLD and XLD manifest different mechanisms that dominate their modes of variability (seasonal vs synoptic). But, how do the synoptic variations in pCO 2sea relate to winddriven upper-ocean physical variability? We address this question first via qualitative analysis of prominent wind events in the time series. During intense wind events (>20 m s −1 , Supplementary  Fig. 2), when the XLD > MLD (particularly during two storm events centred on 28 December 2018 and 5 January 2019, Fig. 2), the mixing vertically entrained the underlying UCDW with its high salinity and DIC into the mixed layer (since both salinity and DIC increase rapidly with depth at the top of the UCDW; Fig. 1c and Fig. 2c). As expected, the entrainment events are followed by increased salinity in the mixed layer ( Fig. 2c), elevated pCO 2-DIC ′ (Fig. 2d), and a reversal from ingassing to outgassing (Fig. 2e) that is indicative of higher DIC in the mixed layer. These two events are associated with two of the largest CO 2 outgassing peaks observed during the experiment (Fig. 2e) and shifted the mean sea-to-air F CO2 by 50% from −0.12 mol C m −2 yr −1 (i.e. the mean excluding these two events of positive F CO2 ) to −0.06 mol C m −2 yr −1 (i.e. the mean of the full record in Fig. 2e including these two events). The two mixing events described above suggest that stormdriven mixing and entrainment provide a significant driver of transient outgassing. However, consideration of a later storm event reveals a more nuanced perspective. During the strong wind event on the 28 of January 2019, when the XLD > MLD (Fig. 2b) and the friction velocity reached its highest value during the multiglider deployment (Fig. 2a), there was only a small increase in mixed-layer salinity (Fig. 2c) and pCO 2-DIC ′ (Fig. 2d). However, the MLD had shoaled and was decoupled from the more saline and DIC-rich UCDW. This suggests that sub-seasonal entrainment of UCDW is not only dependent on the XLD exceeding the MLD but rather on the XLD exceeding the winter MLD maximum (MLD max , located at the 27.3 kg m −3 isopycnal, Fig. 2b), which ultimately sets the depth of the subsurface reservoir of DIC and salinity until the following winter. Thus, the associated wind forcing by passing storms after winter and into late summer requires sufficient intensity to reach this subsurface reservoir of DIC and result in CO 2 outgassing events. It is therefore plausible that progressively stronger buoyancy forcing in summer (e.g.  Fig. 1f. Finally, we observed several other synoptic-scale reversals in ΔpCO 2 between uptake and outgassing events that were not explained by enhanced wind, mixing and entrainment. Although the ΔpCO 2 gradient was reversed during such events, the magnitude of the CO 2 outgassing flux was considerably less than the two entrainment events explained above due to the low wind speed and small K w in Eq. (1) (Fig. 2e). The question that therefore arises is what could explain these further synoptic variations in pCO 2-DIC that do not obviously coincide with strong wind and entrainment? The following section provides insight to this question by invoking a dynamical model for wind-driven synoptic variability of pCO 2-DIC ′ that is driven by both vertical entrainment and meridional Ekman transport.
A dynamical model for synoptic variability of surface ocean pCO 2 . We establish a model that combines two upper-ocean responses associated with the passage of a storm: storm-driven vertical entrainment, which is defined here to occur when the XLD exceeds the MLD max, and storm-driven meridional advection by the Ekman flow (Ekman advection), which drives the meridional displacement of waters in the wind-driven mixing layer. Before proceeding to the physical model, it is important to note that pCO 2sea is not a conservative tracer and it can be estimated from DIC and total alkalinity (A T ) using an empirical function, G 35 . The conceptual framework is to describe synoptic variability of pCO 2-DIC ′ in terms of synoptic variability of conservative DIC and A T : Y Ek (DIC, A T ) represents the anomalies of DIC and A T due to north-south displacements of water masses from the timeintegrated Ekman advection: Where the meridional Ekman velocity is defined by v Ek ¼ À τ x ρ sw XLDf , τ x is the zonal surface wind stress, ρ sw is the reference density of seawater, and f is the Coriolis parameter (or inertial frequency which is negative in the Southern Hemisphere by convention), XLD is the mixing-layer depth, which is assumed equal to the Ekman depth.
is approximated by the meridional gradient of DIC and A T computed using climatological concentrations of A T and DIC derived empirically 36,37 (dashed line, Fig. 3a, b). The temporal variability of the gradients is neglected for simplicity, based on the modest spread in the observed gradients at a few different times (Fig. 3b) and the strong performance of the model (Fig. 3d). The synoptic anomalies of Y Ek were separated by removing the 10-day rolling mean (Supplementary Fig. 7).
Z ENT (DIC, A T ) represents the DIC and A T anomalies due to time-integrated vertical entrainment, which is irreversible and non-negative (unlike Ekman advection). That is, where using and MLD max and XLD are positive by convention. It may be noted that this model omits biological sources and sinks as well as many transport processes, including diffusion through MLD, lateral mixing and several lateral and vertical advective processes, all of which turn out to be less significant than Ekman advection in driving the observed synoptic variability of pCO 2-DIC . In Eq. (4), MLD max is the maximum MLD, C surf and C deep are averaged concentrations of DIC or A T within the surface to MLD max and in the subsurface between the MLD max to~20 m below, respectively (see ref. 38 for a derivation of the entrainment tendency term in a surface-layer average heat budget). The time series of DIC and A T anomalies associated with entrainment are calculated numerically and added to the Ekman anomalies to obtain the surface ocean pCO 2 anomaly (see Methods for more details on the model).
To anticipate the model behaviour, consider the input horizontal and vertical DIC and A T gradients in Eqs. 3-4 ( Fig. 3a, b; Fig. 1c, d), in addition to the wind forcing and mixing reported above (Fig. 2, Supplementary Figs. 5 and 6). In particular, the subpolar Southern Ocean has vertical and meridional gradients in DIC and A T , with higher DIC and A T at depth and further south (Fig. 1c, d). The increasing surface concentrations of DIC and A T to the south (Fig. 3b) are primarily a result of the large-scale meridional upwelling of UCDW (high in DIC and A T ) in the south (Fig. 1c, d) and the effect of southward decreasing sea surface temperature on CO 2 solubility 39 . At 54°S, 0°E, the DIC lateral and vertical gradients were stronger than A T (Fig. 3a, b). Thus, the storm-driven supply of DIC into the mixed layer is more sensitive than A T , which is important given that DIC and A T impact pCO 2sea in opposing ways. This is illustrated by Fig. 3c, a DIC-A T vector plot 40 , which compares the impact of Ekman and entrainment transports on the relative change of DIC and A T . The ratio of A T /DIC sets the sensitivity of the pCO 2sea to both Ekman and entrainment processes (Fig. 3c). If the ratio of A T /DIC were about 1:1, the slope would be along the pCO 2sea isolines and the impact of Ekman advection and vertical entrainment on pCO 2sea would be negligible (Fig. 3c). Our model shows that the A T versus DIC slope (m) for entrainment is 0.3 and for Ekman transport is 0.6, revealing that the pCO 2sea anomalies are positively correlated with DIC anomalies for both entrainment and Ekman processes, but more sensitive to DIC anomalies due to entrainment than to Ekman (Fig. 3c and Supplementary Fig. 8). Thus, the horizontal and vertical gradients of DIC and A T are important factors in determining how efficient these physical processes are in driving the synoptic variability of pCO 2sea .
The resulting time series of the estimated pCO 2-DIC ′ computed using the model (Eq. 2) is compared to the Wave Glider observed pCO 2-DIC ′ (Fig. 3d). Strikingly, the modelled pCO 2-DIC ′ reproduces most of the observed synoptic variability in pCO 2-DIC ′, accounting for 60% of the total observed variance in pCO 2-DIC ′. However, there are some discrepancies between the estimated and the observed values, such as the phasing of the estimated variability is not always aligned with the observed pCO 2-DIC ′ (e.g. events on 25/12 and 05/02 had slightly delayed observed responses). Likewise, the estimated magnitude of the response is sometimes underestimated (12/01) or overestimated (30/12).
Nevertheless, we conclude that physical transport associated with Ekman advection and entrainment and encapsulated in the model is the dominant cause of the observed synoptic variability in pCO 2-DIC , ΔpCO 2 and pCO 2sea . All other transport and biological processes likely explain less than 40% of the variance and are therefore subdominant. For example, variability in biological sources and sinks of DIC may explain a small fraction of the synoptic variance in pCO 2-DIC , but the amplitude of synoptic variations in net primary productivity derived from in situ optical measurements are estimated to be an order of magnitude too weak to explain the observed synoptic variability in pCO 2-DIC ( Supplementary Fig. 9).
The model also allows for the separation of Ekman advection (grey shading) and entrainment (blue shading) contributions to pCO 2-DIC ′ variability (Fig. 3d). We find that storm-driven Ekman displacement dominates entrainment and explains most of the synoptic pCO 2-DIC , ΔpCO 2 and pCO 2sea variations during the deployment, including during the strong entrainment events discussed previously (Fig. 3d). But entrainment does cause rare large pCO 2-DIC ′ anomalies and it contributes substantially to the strongest outgassing fluxes F CO2 observed on 28 December 2018 and 5 January 2019 that result from a synergistic combination of positive pCO 2-DIC ′ and ΔpCO 2 due to Ekman advection and entrainment as well as large K w from strong winds (Fig. 2d).

Discussion
In this study, we present observations of the coupled physicalcarbon processes by which storms drive synoptic variability of ΔpCO 2 and brief, but strong, outgassing events in the subpolar SE Atlantic. Hence, the results from this process study raise important broader questions: how prevalent is this strong synoptic variability in the coupled ocean physics-carbon system around the entire subpolar Southern Ocean? And what does it mean for the larger spatial and seasonal-inter-annual CO 2 flux dynamics?
Prevalence of storm-driven synoptic variability in pCO 2 across the subpolar Southern Ocean. In order to shed light on the circumpolar prevalence of regions of synoptic variability and the larger-scale implications as well as to test the proposed Entrainment-Ekman conceptual model for pCO 2-DIC , we have applied the model (Eqs. 2-4) over the dynamically comparable circumpolar upwelling zone of the Southern Ocean (between the Polar Front and the northern limit of sea ice in winter; Fig. 4; see Methods). The magnitude of the estimated synoptic variance (i.e. mean square anomaly) of pCO 2-DIC can reach 25 μatm 2 in places with a mean-variance of 6 μatm 2 (Fig. 4a). For perspective, this 6 μatm 2 variance due to storms is approximately half in the magnitude of the summer mean inter-annual variance of pCO 2sea estimated for the subpolar region of 13 μatm 2 computed from the CSIR-ML6 observation-based product 1,41 . Moreover,~30% of the surface area of the subpolar ocean has a synoptic variance that is greater or equal in magnitude to the summer mean inter-annual variance of pCO 2sea , thus synoptic variability of pCO 2sea is potentially a widespread dominant mode of variability in the Southern Ocean.
To assess the robustness of the estimated spatial distribution of synoptic variability of pCO 2-DIC ′, we spatio-temporally collocate the estimate derived from Eqs. 2-4 with recent Saildrone Antarctic circumpolar pCO 2 observations 28,42 . We find a statistically significant agreement between our estimate of the pCO 2-DIC ′ variance and the observed spatio-temporal synoptic variability of pCO 2-DIC ′ (r 2 = 0.4, Fig. 4a) from the Saildrone dataset. Regions of higher variance in the modelled pCO 2-DIC ′, such as in part of the SE Atlantic, coincide with higher observed pCO 2-DIC ′ variance observed by the Saildrone (indicated by the size of the hexagon). In the western Pacific sector, low estimated synoptic-frequency variance overlapped with lower measured variance (Fig. 4).
In agreement with the in situ observations from the multiglider deployment (see Fig. 3), lateral Ekman advection over vertical entrainment is the dominant physical driver of the estimated high-frequency variability explaining about 92% of the model estimated variance of pCO 2-DIC ′ across the subpolar Southern Ocean (Fig. 4b). Importantly, the spatial variation in the meridional gradients of pCO 2-DIC explained most of the spatial variation in high-frequency temporal variability of pCO 2-DIC ′ (Fig. 4a, c, d). This is consistent with Ekman advection scaling with the zonal wind stress and meridional gradients of DIC and Total Alkalinity ∂ DIC; A T ð Þ ∂y , because the spatial variability of the high-frequency wind is relatively uniform across this region (synoptic atmospheric variability occurs across large spatial scales of order 1000 km, refer Fig. 4c).
The hypothesized dominance of Ekman advection of mean gradients as a driver of the synoptic variability of pCO 2-DIC has some somewhat surprising and important implications. First, it implies that time-mean gradients on the timescales relevant to the synoptic lateral advection (that is days to weeks). This result in turn implies that the large-scale (>500 km) A T and DIC fronts are fairly stable in time (e.g. as seen by the low meridional variability of the DIC and A T observed in Fig. 3b). Drivers of large-scale variability such as changes in the large-scale circulation, biological productivity and air-sea fluxes evidently do not cause substantial inter-annual or even seasonal deviations from the time-mean . Stirring by mesoscale eddies is relatively ineffective at producing A T and DIC anomalies (e.g. via frontogenesis, see ref. 43 ) compared to the mechanisms that are destroying the mesoscale gradients. Relatedly, the dominance of Ekman advection also implies that Ekman velocities dominate all other sources of synoptic velocity variability in the mixed layer, that is mesoscale and submesoscale turbulent velocities are relatively weak compared to Ekman velocities on synoptic timescales. A full explanation of these results and a broader evaluation of these hypotheses is beyond the scope of this work, but it is important to recognize that the explanatory power of Ekman advection of mean gradients and hence the accuracy of the extrapolations in Fig. 4 depend on the relative weakness of both variability in meridional gradients of pCO 2-DIC and non-Ekman synoptic velocities, as inferred from observations in the SE Atlantic multiglider deployment. We have not investigated what sets the magnitude of these large-scale mean lateral gradients of pCO 2-DIC (Fig. 4d) or the mesoscale kinetic energy; however, these topics have been explored in other papers and are thought to be associated with the large-scale ocean fronts of the Antarctic Circumpolar Current (e.g. refs. [44][45][46] ).
Implications of storm-driven synoptic variability for carbon dynamics on longer timescales. Another outstanding question is whether or not the storm-driven Ekman advection and entrainment have implications for the slower seasonal or inter-annual carbon dynamics of the Southern Ocean? Here, we address this question in two parts, first focusing on the Ekman advection and then entrainment. Perhaps the most striking result derived from the glider deployment and the subsequent extrapolation across the subpolar Southern Ocean is that a simple Ekman advection of mean meridional gradients explains the majority of the synoptic variance in pCO 2sea and ΔpCO 2 (Fig. 3 and Fig. 4). However, to a first approximation (and by definition in Eq. 3), oscillatory synoptic perturbations to Ekman advection are reversible and do not sum to impact the mixed-layer DIC budget and F CO2 on longer timescales. To this level of approximation, perturbations to Ekman advection only locally modulate the sign of ΔpCO 2 where the seasonal drivers associated with thermal and non-thermal effects yield small mean |ΔpCO 2 | relative to the synoptic perturbations, and thus the sign of ΔpCO 2 and the resulting F CO2 is highly sensitive to the synoptic perturbations and varies on synoptic timescales. In addition, the observed mean ΔpCO 2 is only sensitive to the synoptic perturbations in Ekman advection to the degree that the synoptic variance is also relatively large compared to the number of synoptic events or the duration of sampling (i.e. standard errors are large). However, oscillatory Ekman advection can be rectified in other subtle ways that are not captured at this level of approximation. For example, oscillatory advection of otherwise static ocean pCO 2sea spatial gradients under a spatially and temporally variable atmosphere may modify the average F CO2 over longer time intervals due to correlation between Ekman-driven ΔpCO 2 anomalies and wind speed as well as the non-linear dependence of F CO2 on wind speed. In addition, the combination of oscillatory Ekman advection and intermittent ocean mixing (e.g. during entrainment events) may induce lateral mixing via shear dispersion 47 that irreversibly sums to impact the slower evolution of the mixed-layer pCO 2sea . However, it is beyond the scope of this manuscript to quantify these rectified effects of Ekman advection for the large-scale dynamics of DIC and the air-sea CO 2 flux.
On the other hand, entrainment events, which have been shown to have a much smaller contribution to the synoptic variance than Ekman advection (Figs. 3 and 4a, b), connect the subsurface carbon-rich UCDW to the mixed layer irreversibly. Thus, all entrainment events sum to impact the mixed-layer DIC budget and hence F CO2 on longer timescales (seasonal to interannual). The probability of sampling short storm-driven entrainment events like those observed by the multi-glider deployment with a 10-day (e.g. floats) or greater sampling period (e.g. ships) is very low and the response to entrainment is obscured by Ekman advection in any case. Thus, it is difficult to observationally quantify intermittent synoptic entrainment fluxes across the entire subpolar Southern Ocean as we do in the SE Atlantic with this data from paired gliders; coarser spatiotemporal sampling may alias this variability 17,28 . Hence, we use the model Eqs. 2-4 to provide an estimate of the magnitude of the time-averaged synoptic entrainment flux across the subpolar Southern Ocean. Figure 5 quantifies the annual mean entrainment flux of DIC (Eq. 4, see also Methods) and compares it (for perspective) with the magnitude of the climatological seasonal amplitude of F CO2 (see also Fig. 1a, b). It shows that mean storm- driven entrainment flux (~3.5 mol C m −2 y −1 ) is of a similar order of magnitude to the amplitude of the seasonal cycle in F CO2 (~2.1 mol C m −2 y −1 , peak to trough) as well as the time-mean F CO2 (Fig. 1a;~−0.1 mol C m −2 y −1 ). Hence, even small variations in the synoptic entrainment flux of DIC have a large impact on the mixed-layer DIC budget relative to F CO2 . If changes in synoptic entrainment go uncompensated by changes in other sources/sinks of mixed-layer DIC such as biological export production (which is plausibly of the same order of magnitude and opposite sign; see Supplementary Fig. 4) or other physical transport processes, the entrainment will drive changes in DIC and F CO2 . Consideration of the spatial structure of synoptic entrainment in Fig. 5 shows that entrainment exhibits substantial spatial variability and is particularly strong in the South Atlantic where the MLD max is comparatively shallow relative to the Pacific basin (ref. 22 their Fig 10c) and storms are more frequent and stronger 13 . The spatial variability of the synoptic entrainment highlights the variable circumpolar implications of the observation reported above that storm-driven entrainment is sensitive to the winter MLD maximum (which sets the depth of UCDW reservoir) (Figs. 1b, c and 2b) and to stabilizing buoyancy forcing (which prevents storm-driven vertical mixing from reaching the UCDW reservoir during later summer months, e.g. in Fig. 2). Finally, these results emphasize that entrainment, which depends on the vertical gradients of DIC and A T and the MLD (Eq. 4), exhibits a quite different spatial structure than synoptic variance due to Ekman advection, which depends on the meridional gradient in pCO 2-DIC (Eq. 3). Nevertheless, in both processes, different underlying oceanic conditions (in addition to the atmospheric conditions) are crucial determinants of the oceanic response to storms. Although a complete analysis of the mixed-layer DIC budget is beyond the scope of this paper, the results in Fig. 5 indicate that it is possible that storm-driven entrainment cumulatively impacts the mean pCO 2sea and CO 2 flux and influences inter-annual 28 and spatial variability of CO 2 outgassing ( Fig. 1 and Fig. 5) through the interactions between annual changes in stormcharacteristics 48 , the seasonal cycle of the mixed layer, and variations to the depth of the winter MLD maximum 49 . We leave tests of these hypotheses to future work. The success of the conceptual model (Eqs. [2][3][4] in describing observations in the Atlantic sector (Fig. 3) coupled with the spatially variable and significant implications of storms for the mixed-layer DIC budget on a range of timescales (Figs. 4 and 5) motivate future field experiments and comprehensive modelling that can accurately quantify and predict the physical-carbon dynamics of the ocean mixed layer down to synoptic timescales more broadly around the subpolar Southern Ocean.

Methods
Experimental design and region of interest. The basis for the observations conducted in this study is the SOSCEx-STORM experiment, which aims to simultaneously observe the passage of storms, their associated surface wind stress and corresponding upper-ocean response in physics and CO 2 . The SOSCEx-STORM experiment forms part of a larger observational program, the Southern Ocean Seasonal Cycle Experiment (SOSCEx), detailed in ref. 50 . In SOSCEx-STORM, twinned gliders, a surface Wave Glider and a profiling glider, were programmed to sample together, allowing for a high-resolution view of the coupled atmosphere and upper-ocean processes. The field site chosen was at 54°S, 0°E, which is located south of the Antarctic Polar Front and in the globally significant outgassing sector of the Southern Ocean (Fig. 1). The platforms sampled for 2 months (56 days), between 18 December 2018 and 12 February 2019.

Autonomous observing platforms
Wave Glider integrated with surface CO 2 sensor. The Liquid Robotics SV3 Wave Glider (WG) was fitted with an Airmar WX-200 Ultrasonic Weather Station mounted on a mast at 0.7 m above sea level and sampled wind speed and direction at a rate of 1 Hz, averaged into 10 min bins. The surface winds were corrected to a height of 10 m above sea level as in ref. 51 . The WG was equipped with a SeaBird Glider Payload CT-cell, measuring surface ocean temperature and conductivity at 1 Hz, averaged into 20 min bins. In addition, the WG was fitted with a VeGAS-pCO2 (Versatile Glider, Atmospheric and Ship pCO 2 high Precision pCO 2 analyzers) measuring atmospheric and ocean pCO 2 . The VeGAS-pCO 2 sensor is based on the well-established NDIR (Licor -Li-820) linked equilibrator units 52,53 but with a significant redesign to improve accuracy (<1 µatm), precision (<1 µatm) through more effective drying and temperature control, equilibrator design and long term stability that also reduced the frequency of reference gas calibration from every sample to every 2 h. The unit was installed and linked to the SV3-WG control unit which enables remote communication and to send real-time data. These instruments have just recently been successfully assessed in the ICOS Ocean Thematic Centre instrument intercomparison study and those results will be published through ICOS. Outlier data points were removed from the WG pCO 2 by applying a global cut-off upper 99.9 percentile and lower 0.1 percentile of the discrete temporal difference. A rolling mean with a window size of half the local inertial period (i.e. about 8 h) was applied to pCO 2 WG data.
Slocum integrated with an RSI MicroRider. The water column observations were collected using a Teledyne Webb Slocum G2 glider. This profiling glider was equipped with a pumped SBE conductivity and temperature sensor, and a Rockland Scientific (RSI, Canada) MicroRider for microstructure measurements. The MicroRider was fitted with two piezo accelerometers and two air-foil shear probes and only collected data during the Slocum climbs. This choice was motivated by the need to increase the battery endurance by sampling for only half of the duration (about 3 months), and the choice of climbs over dives was made to ensure dissipation estimates as close to the surface as possible. The Slocum dataset was processed with the GEOMAR MATLAB toolbox and post-processing gridding (i.e. vertical interpolation to 1 m bins) and quality control of the glider data was carried out with GliderTools 54 . A Savitzky-Golay 55 filter was applied to the glider salinity and temperature profiles to remove spurious spiking in the data (as recommended in GliderTools) and smoothed further with a rolling mean (9 profiles or about 18 h window). The GEOMAR toolbox includes a hydrodynamic glider flight model that produces a time series of flow past the sensor and angle of attack (AOA) which are required for the processing of accurate microstructure measurements (ref. 56 , section 4). The microscale velocity shear was obtained from the shear probes measuring the two orthogonal components of the shear along the axis of the instrument ∂v ∂x À and ∂w ∂x Á . The MicroRider accelerometers, which obtain high accuracy measurements of mechanical-driven or impact-driven vibrations felt by the Slocum glider, were used to remove vehicle vibration contamination of the shear data 57 . The shear data from both probes were then processed into dissipation rate (ε) estimates using a four second FFT length and 12 second averaging 56 . The processing of the microstructure data was based on the routines provided by RSI (ODAS v4.04 software). The dissipation rate values were calculated assuming isotropic turbulence ε ¼ 7:5υ ∂w ∂x À Á 2 (here written for the w component) where υ is the seawater viscosity. The small-scale shear variance (the term with the overline) was obtained by integrating the wavenumber spectrum of shear in a wavenumber range that is relatively unaffected by noise and corrected for the unresolved variance using the empirical model from ref. 58 . The dissipation estimates underwent further quality controls. This included, globally, any segments with AOA larger than 6°, pitch larger than 30°, and Figure of Merit (a measure of spectra fit the Nasmyth spectrum, high values are a poorer fit) larger than 2 are excluded 56 (RSI Technical note 039, https://rocklandscientific.com/support/knowledge-base/technical-notes). We further excluded individual dissipation estimates which seemed unrealistic in value (e.g., data spikes due to plankton interactions on the shear probes).

Derived metrics
Friction velocity and theoretical dissipation. Wind observations from the Wave Glider and dissipation estimated from shear probes on the MicroRider are used to investigate the link between wind and ocean turbulence via the friction velocity u * defined as, u Ã ¼ , where ρ sw is the density of seawater and τ is the wind stress estimated using 59 . Similarity theory states that wind, through u * , impacts the dissipation rate of turbulent kinetic energy (ε) via ε ¼ uÃ 3 kz , where k = 0.41 is the von Karman's constant and z is the depth range below the surface ocean 60 .
Mixed-layer depth and Mixing-layer depth. The mixed-layer depth, MLD, was calculated as the depth from the surface where the density first exceeds its surface value by 0.03 kg m −3 as in refs. 61,62 , using the glider observed density sections. The mixing-layer depth, XLD, was estimated using the MicroRider derived ε profiles, as the depth from the surface where ε first drops below 10 −8 m −2 s −3 , following 33 .
Wave Glider estimated A T , DIC and CO 2 flux. The WG observed surface ocean temperature and salinity were used to estimate A T 63 , DIC was derived from estimated A T and WG observed pCO 2sea using PyCO2SYS 64 (with options for equilibrium constants as follows: K1 and K2 are from Mehrbach, refit by Dickson and Millero, KSO4 from Dickson and TB from Uppstrom). Uncertainties in derived A T and DIC are about 5-10 µmol kg −1 . WG observed temperature, salinity and wind speeds along with the Japanese 55-year Reanalysis (JRA-55-do 65 ) atmospheric sea level pressure (Supplementary Fig. 2) were used to compute air-sea CO 2 fluxes using the bulk formulation with python package Seaflux.1.3.1 (https://github.com/ lukegre/SeaFlux) 66 .
Net primary productivity. Net primary productivity (NPP) is estimated from satellite ocean color (Ocean Colour-CCI 67,68 ) using three different primary models: the Carbon-based Productivity Model (CbPM 69 ), the Vertically Generalized Production Model (VGPM 70 ) and the Platt model 71 . In addition, NPP is estimated from bio-optical measurements on the Slocum glider, which was fitted with a WETLabs ECO puck™ (BB2Fl-470/700), thus measuring chlorophyll-a fluorescence (proxy for phytoplankton concentration) and two wavelengths of optical backscattering by particles, bbp(470) and bbp(700). The optics data were cleaned and processed following procedures recommended by GliderTools 54 . Backscatter, bbp(700), was converted to phytoplankton carbon following Behrenfeld et al. 72 . Slocum glider NPP was estimated from the phytoplankton carbon, chlorophyll and collocated photosynthetically available radiation (PAR) from MODIS (the Slocum glider did not have a PAR sensor) following the CbPM 69 using python code Pri-maryProductionTools.py (https://github.com/isgiddy/roammiz-seaice-impactsorganic-carbon/blob/v0/src/PrimaryProductionTools.py) 73 adapted from Arteaga et al. 74 .
Generalization of the conceptual model across the subpolar outgassing region. We have applied the box model for high-frequency variability of surface ocean pCO 2 (defined in Eqs. 2-4) across a dynamically representative zonal band, which we define to be constrained to (1) a region that is not impacted by sea ice, (2) a region of mean outgassing, (3) a region with similar mean wind forcing and (4) mean mixed-layer depths. We used 3-hourly reanalysis winds to compute u * using winds from the Japanese 55-year Reanalysis (JRA-55-do 65 ); monthly mean MLD were estimated using a density threshold of 0.03 kg m −3 from density field derived from EN4.2.1 interpolated fields of temperature and salinity 75 ; and temporal varying (i.e. monthly mean) lateral and vertical gradients from the climatology of A T and DIC were taken from refs. 36,37 . We translate wind variability via friction velocity (u * ) into mixing-layer depth (XLD) variability assuming the strong relationship between u * and XLD holds true for these dynamically comparable regions ( Supplementary Figs. 5 and 6). Intermittent high-frequency supplies of DIC and A T due to high-frequency Ekman advection and wind-driven vertical entrainment were computed following (Eqs. [2][3][4] and the anomalies of DIC and A T were iteratively added to baseline climatological means of DIC and A T . Changes in temperature were not included as the focus was on non-thermal drivers. Finally, PyCO2SYS 64 was used to compute pCO 2-DIC from the physically-driven changes in A T and DIC to generate high-frequency temporal variability of pCO 2-DIC anomalies across the Subpolar Southern Ocean (as shown in Fig. 4).
It is important to recognize that the generalization of the conceptual model (Eqs. 2-4) across the subpolar Southern Ocean for the full year ( Fig. 4 and Fig. 5) is not validated by in situ observations at locations and times beyond those isolated validations reported here, and there are several assumptions that underpin the model. In particular, the model omits many potentially significant processes to isolate those that we find to be most important in the SE Atlantic. For example, non-thermal (pCO 2-DIC ) components of pCO 2sea are also driven by processes that include net community production (NCP), calcification, air-sea CO 2 exchange and freshwater fluxes, as well as various advective and mixing processes other than Ekman and entrainment, the influence of which is not included in the model. The premise for these omissions is that these terms have a comparably smaller impact on pCO 2-DIC on these shorter 1-10-day timescales, as observed in the SE Atlantic during summer (Fig. 3d, Supplementary Fig. 9). Particular caution must be exercised when using the results of the conceptual model in regions where the neglected factors may be more significant, either because the strength of one or more of the neglected factors is greater or because the variability driven by Ekman and entrainment is weaker. With regard to the neglected NCP, for example, NPP is much stronger and more variable than at 54°S, 0°W in some localized regions surrounding subpolar islands and off the coast of South America ( Supplementary  Fig. 4). Conversely, physical DIC variability due to Ekman advection is substantially weaker in areas where meridional DIC gradients are weak (Fig. 4d). In addition, even the included entrainment process may not be well represented by the model at all times. For example, neglected effects of convection may have a stronger influence during winter than is reflected in the empirical relationships based on summertime observations ( Supplementary Figs. 5 and 6). Furthermore, the magnitude and vertical extent of upper-ocean mixing and entrainment are influenced by other physical processes such as wave action 76,77 and submesoscale dynamics [78][79][80] , which are not included. Finally, the generalization of our conceptual model is dependent on the accuracy of both reanalysis winds (JRA-55do reanalysis has been shown to perform well in this region with a mean difference of −0.02 m s −1 ± 0.8 m s −1 with the WG wind speeds and refer to Supplementary  Fig. 2 65 ) and on statistically inferred products of climatological DIC and A T 36,37 .