Air‐Sea CO2 Fluxes Localized by Topography in a Southern Ocean Channel

Air‐sea exchange of carbon dioxide (CO2) in the Southern Ocean plays an important role in the global carbon budget. Previous studies have suggested that flow around topographic features of the Southern Ocean enhances the upward supply of carbon from the deep to the surface, influencing air‐sea CO2 exchange. Here, we investigate the role of seafloor topography on the transport of carbon and associated air‐sea CO2 flux in an idealized channel model. We find elevated CO2 outgassing upstream of a seafloor ridge, driven by anomalous advection of dissolved inorganic carbon. Argo‐like Lagrangian particles in our channel model sample heterogeneously in the vicinity of the seafloor ridge, which could impact float‐based estimates of CO2 flux.

transport of anthropogenic carbon (Ito et al., 2010), and that intensified residual upwelling downstream of regional topographic features provides an important conduit for deep, natural carbon to enter the Southern Ocean surface (Brady et al., 2021).Despite the potentially important role that these regional topographic features play in the global carbon budget, no study has directly quantified the influence of seafloor topography on Southern Ocean air-sea CO 2 flux nor addressed the potential effects these features may have on Lagrangian observations of the Southern Ocean.
Here, we use an idealized, high-resolution ocean general circulation and biogeochemical model to assess the role of seafloor topography in Southern Ocean air-sea CO 2 fluxes and the ability to quantify these fluxes via Lagrangian observations.Our study demonstrates that seafloor topography has a substantial impact on local CO 2 flux via topography-driven advection of dissolved inorganic carbon (DIC).Lagrangian particles tend to heterogeneously sample the surface pCO 2 in the vicinity of topography, and this can affect estimates of average air-sea CO 2 fluxes over the region.In Section 2, we present the methods used, in Section 3 we present the results.In Section 4 we discuss and conclude.

Model Description
For this study, we use an idealized-geometry MITgcm ocean channel model (Youngs & Flierl, 2023) and couple it to a simple ocean biogeochemical model (Dutkiewicz et al., 2005;Lauderdale et al., 2016).The channel is 4,000 km long and 2,000 km wide with 10 km horizontal resolution (Figure 1) with a total depth of 4,000 m with 32 vertical levels, from 10 m vertical grid spacing at the surface to 280 m at the bottom.We represent seafloor topography using a 2,000 m tall Gaussian ridge with a 200 km half-width, centered 800 km downstream of the channel entrance spanning the channel north to south (Figure 1).The domain is periodic with the outflow in the east reentering in the western boundary and free-slip walls at the north and the south.The model is integrated using a 600 s time step, an exponentially varying diffusivity (0.01 m 2 s −1 decreasing downward to 1 × 10 −5 m 2 s −1 ), and linear bottom drag with a drag coefficient of 1.1 × 10 −3 m s −1 .The wind stress is a cosine profile with a maximum value of 0.15 N m −2 at the center of the domain and zero wind stress at the sides (Figure S1 in Supporting Information S1).The salinity is set at 35 PSU and not allowed to vary.
We employ the DIC package from MITgcm to represent biogeochemistry in our model (Dutkiewicz et al., 2005;Lauderdale et al., 2016).This model package carries alkalinity, DIC, dissolved organic phosphate, and phosphate as biogeochemical tracers, and represents biological uptake as a function of phosphate and light availability.Phosphate is fluxed vertically with remineralization and sinking (see more in the SI).The calcium carbonate formation is proportional to the organic phosphorous produced in the surface waters following the parameterization of Yamanaka and Tajika (1996), with sinking and dissolution (Dutkiewicz et al., 2005).
The rate of change of carbon in our model can be described by the following equation (Lauderdale et al., 2016) where C T is the concentration of total dissolved organic carbon, κ is the eddy diffusivity tensor,    ∶ is the biological transformation between carbon and phosphorous and   2 is the air-sea CO 2 flux, h is the mixed layer depth, S bio represents the sources and sinks of biogenic soft tissue, and   3 represents the sources and sinks of biogenic carbonate.Note that this equation neglects the dilution by freshwater fluxes, which in our case is appropriate due to a lack of salinity or freshwater forcing.
The model is initialized with a uniform surface ocean pCO 2 of 270 ppm with DIC and alkalinity at the northern boundary sponge region relaxed to prescribed preindustrial DIC and alkalinity profiles based on GLODAPv2.2016(Key et al., 2015;Lauvset et al., 2016) (Figure S3 in Supporting Information S1), and spun up for 30 years for the biogeochemical and physical tracers to reach an approximate steady-state (Table S1 in Supporting Information S1).At the end of the spin-up period, our model simulates similar Southern Ocean-integrated pre-industrial air-sea CO 2 fluxes (0.1 mol m −2 yr −1 ) as those estimated from more realistic model configurations (0.13 mol m −2 yr −1 ) (e.g., Lovenduski et al., 2007).The atmospheric CO 2 is held at a preindustrial value of 270 ppm and the air sea fluxes are calculated via a bulk formula.

Particle Tracking
We model idealized "Argo" float trajectories to estimate how well a biogeochemical Argo float array can sample the air-sea carbon fluxes as a function of float density.We use the Ocean Parcels package to track idealized Argo floats (https://oceanparcels.org/)(Lange & van Sebille, 2017).We release 800 floats spaced uniformly throughout the model domain.Real Argo floats park at 1,000 m depth for 10 days between profiles, so in our simulations the particles are advected using daily-averaged velocities at 1,000 m; they sample the surface ocean pCO 2 at their position every 10 days.Idealized (and real) Argo floats are advected in a 2D velocity field, which could be convergent/divergent because they cannot follow the vertical flow.Idealized floats are advected for 1 or 3 years.We run four collections of experiments: 10 floats for 1 year, 33 floats for 1 year, 100 floats for 1 year, and 33 floats for 3 years.For each collection, we select 100 instances of that number of floats (from a larger collection of 800 floats) for the time series in question to calculate statistics.We use the randomly subsampled float data to create a climatology using objective mapping (e.g., Figure 3b).From the mapped pCO 2 , we calculate the air-sea carbon fluxes using the same equations used by the model (Wanninkhof, 1992).
Objective mapping is a commonly used and well-justified technique for mapping sparsely sampled data to estimate regional averages (Dong et al., 2008;Friedrich & Oschlies, 2009;Reeve et al., 2016).We create climatologies of these samples using the ordinary kriging method with the PyKrige Python package (https://github.com/GeoStat-Framework/PyKrige/).Here, the various terms for the Gaussian variogram are fit using the data from the selected floats to create the most optimal map.

DIC Budget
We investigate the asymmetry of the carbon properties in the channel model.Both air-sea CO 2 flux and surface DIC concentration exhibit large zonal asymmetry, with enhanced CO 2 outgassing and elevated surface DIC located just upstream of the undersea ridge (Figures 1b and 1c).Away from the influence of topography, our model exhibits moderate outgassing of CO 2 near the southern boundary, with weak uptake in the northern part of the domain (Figure 1b), which together contribute to an average flux of about −0.07 mol C m −2 yr −1 .At the latitudes of the topographic ridge, however, we find sea-air CO 2 fluxes that exceed 7 mol C m −2 yr −1 and outgassing that extends to the northern boundary of the domain, with an average flux of 0.8 mol C m −2 yr −1 .The enhanced carbon flux is located in the region where the barotropic flow turns north as it approaches the ridge (Figures 1b and 1c).This region is characterized by elevated surface DIC concentrations relative to the zonal mean for the domain (Figure 1c).We also show that as the wind stress forcing changes, the pCO 2 flux changes are driven by changes in advection of DIC not other terms like temperature forcing or changes in alkalinity (Figure S5 in Supporting Information S1), highlighting the importance of the advection of DIC.
We investigate the drivers of the elevated surface ocean DIC upstream of the topographic ridge by quantifying the terms in Equation 1 averaged over the top 50 m.DIC advection tends to increase DIC upstream of the ridge, while sea-air CO 2 flux tends to decrease DIC in this same region (Figures 2a and 2b).In contrast, biological productivity tends to decrease DIC relatively uniformly over the domain, with only a slightly larger influence upstream of the ridge, and DIC diffusion exhibits only a small influence on upper ocean DIC tendency across the domain (Figures 2c and 2d).The elevated net DIC advection upstream of the ridge is mostly driven by vertical advection (Figure S4 in Supporting Information S1), though the contribution from the horizontal advection of DIC is non-negligible, especially in the northern portion of the model domain (Figure S4 in Supporting Information S1).Thus, results from our DIC tendency budget suggest that enhanced vertical advection of DIC upstream of the ridge is responsible for the locally elevated DIC, and by inference, the enhanced outgassing of CO 2 in this region.Our model also simulates elevated sea-air CO 2 flux and surface ocean DIC in the northern portion of model domain over the ridge, albeit with lower magnitudes than in the region upstream of the ridge (Figure 1).Here, the elevated DIC is driven by DIC advection (Figure 2), with horizontal DIC advection playing a key role (Figure S4 in Supporting Information S1).

Sampling Heterogeneous Carbon Fluxes
Topography-induced heterogeneity may challenge observation of ocean carbon processes.We quantify the ability of autonomous, Lagrangian floats to sample surface ocean DIC and associated CO 2 fluxes by adding idealized particles to our model domain.We test four deployment strategies (a) 10 floats for one year, (b) 33 floats for one year, (c) 100 floats for one year, and (d) 33 floats for 3 years.For each float number and duration we select 100 collections of random initial conditions.We calculate the error by subtracting the model truth from the calculated air-sea CO 2 fluxes, integrating over the residual and normalizing by the integrated value of the model truth air-sea CO 2 fluxes.As such, our error estimate is fairly conservative; the error would certainly be larger using a square error metric.
Our idealized sampling approach reveals substantial biases in the domain-integrated CO 2 flux, as compared to the model truth.With 10 floats, the interquartile range of the air-sea CO 2 flux error is large, from a 113% overestimate to a −146% underestimate, with larger extremes in the upper and lower 25% of the realizations.In this  While previous studies have identified topographic upwelling and outgassing of CO 2 , most have suggested that this effect is due to the high eddy kinetic energy downstream of topographic features (Barthel et al., 2022;Brady et al., 2021;Dufour et al., 2015;Tamsitt et al., 2017;Viglione & Thompson, 2016;Yung et al., 2022).Here, by using an idealized process study, we have identified the standing eddy upstream of topography as the dominant cause of upwelling.Future work should test the sensitivity of this result to more complex topography, more realistic atmospheric forcing, and model resolution.This particular mechanism, driven by a standing eddy projection upstream of the domain, may be somewhat reproducible without resolving mesoscale eddies, however, the eddies are necessary for maintaining appropriate gradients.
Our findings suggest that Lagrangian floats may undersample topographically induced biogeochemical anomalies (e.g., DIC, oxygen, nitrate).Future efforts in observational network design should consider alternate methods to estimate the biogeochemistry of topographically influenced regions.One approach is to use alternative technologies such as gliders (e.g., Dove et al., 2021).This study uses the current standard Gaussian objective mapping technique to map surface ocean pCO 2 and infer air-sea CO 2 fluxes.A complementary approach to confronting the challenges posed by Lagrangian autonomous sampling platforms is developing mapping techniques that account for heterogeneous environments such as techniques that utilize information about correlation length scales (Chamberlain, 2022), and those that use ancillary data such as temperature and salinity to map biogeochemical variables (A.Gray, pers. comm.).Such approaches may improve the sampling error in topographically influenced regions.
The idealized model geometry used in this study has enabled mechanistic insights into the drivers of outgassing hotspots at topographic features in the Southern Ocean (Brady et al., 2021;Tamsitt et al., 2017).The insight that barotropic effects have a primary role in driving outgassing hotspots has direct implications for observing system design.Increasing model complexity through more complex and realistic model geometry, improved realism of multiple biogeochemical tracers, finer resolution model configurations, and seasonal variability that can improve representation of wind-current interactions (Kwak et al., 2021) may enable additional insights about the ways that zonal asymmetry influences the Southern Ocean carbon cycle and the coupling between DIC and other biogeochemical factors in the Southern Ocean.
Seafloor topography induces anomalies in both the flow and the surface ocean DIC concentration, leading to sub-optimal sampling of a key region for Southern Ocean CO 2 flux.Through the mechanistic insight provided by this study, we suggest that the current SOCCOM float array has most likely undersampled (rather than oversampled) potential areas of CO 2 outgassing in the Southern Ocean, which could further amplify the differences in CO 2 fluxes estimated from SOCCOM floats and those estimated from ship-or satellite-based observations (Bushinsky et al., 2019;Gray et al., 2018;Long et al., 2021).Topographically influenced regions in the Southern Ocean should be a focus for future biogeochemical observation and modeling programs.

Figure 1 .
Figure 1.The model is a re-entrant channel forced with both a zonal wind and a relaxation to a meridional temperature gradient.Barotropic streamlines are shown with black contours on the top faces.Shading shows temporally-averaged (a) surface eddy kinetic energy, (b) surface carbon dioxide flux, and (c) dissolved inorganic carbon (DIC) concentration.In (c) the right edge shows temporally and zonally averaged DIC concentration with temperature (density) contoured in black.The model geometry is shown in (c).A 2000 m tall undersea Gaussian ridge is centered at x = 800 km.

Figure 2 .
Figure 2. The drivers of the rate of change of DIC (     ; mmol m −3 yr −1 ), as in Equation 1, averaged over the 20 years simulation and the top 50 m: (a) DIC advection, (b) sea-air flux of CO 2 , (c) biology, and (d) DIC diffusion.The vertical lines indicate the location of the top of the ridge.

Figure 3 .
Figure 3. Percent error in the domain-integrated sea-air CO 2 fluxes with Argo-like model sampling for 10, 33, and 100 floats advected for one year and 33 floats advected for 3 years, respectively.The results of 100 trials with randomly initialized floats are shown for each float density.The boxes show the interquartile range and median (orange line) and the whiskers show 1.5 times the interquartile range over the 100 trials.Positive numbers represent anomalous CO 2 outgassing in the float estimate.

Figure 4 .
Figure 4. Air-sea CO 2 fluxes over one year derived from the idealized channel model.(a) Modeled fluxes, (b) fluxes as sampled by 33 randomly spaced particles, and (c) the difference between the sub-sampled fluxes and the model truth, with a 19% overestimate of the fluxes.Gray contours indicate barotropic streamlines, while black lines show the tracks of the 33 floats used to generate the images in panels (b) and (c).