Localizing the Southern Ocean Biogeochemical Divide

The meridional overturning circulation consists of an upper and lower cell. The Southern Ocean Biogeochemical Divide (SOBD) is the boundary between the two cells in the surface of the Southern Ocean, but its location is poorly constrained. Localizing the SOBD is important because biological nutrient utilization north and south of the SOBD have fundamentally different consequences for global ocean primary production and carbon sequestration. Here, we aim to localize the SOBD by releasing virtual Lagrangian particles south of 40°S in an eddying ocean sea‐ice model and compare simulation results with observations. We find that the SOBD is a circumpolar band, where different sectors are shaped by different oceanographic features: (a) Ekman transport, (b) the γ = 27.6 kg m−3 neutral density outcrop, and (c) fronts associated with the Antarctic Circumpolar Current. Our findings help to understand how nutrient utilization in different parts of the Southern Ocean affects the biologically driven carbon sequestration.


10.1029/2022GL098260
3 of 9 mostly south of the SOBD. Thus, a large fraction of the nutrients are trapped in the lower cell and do not fuel primary production elsewhere in the global oceans.
Like the nutrients in the lower cell, nutrients in the upper cell are only partially consumed by phytoplankton during their residence in the surface ocean (i.e., before being subducted within AAIW/SAMW). In contrast to the lower cell, however, the preformed nutrients entrained within AAIW/SAMW remain relatively shallow and eventually become upwelled to fuel primary production and associated atmospheric CO 2 drawdown outside the Southern Ocean (Hauck et al., 2018;Sarmiento et al., 2004). Thus, upper cell primary production in the Southern Ocean competes with downstream primary production outside the Southern Ocean for nutrients. For example, if more upper cell nutrients are utilized, downstream production will be reduced Marinov et al., 2006;Primeau et al., 2013). Indeed, model simulations suggest that as much as 75% of marine primary production north of 30°S is fueled by nutrients supplied by AAIW and SAMW (Sarmiento et al., 2004).
The relevance of the SOBD for the marine nutrient and carbon cycle has been assessed by biogeochemical modeling (Holzer et al., 2014;Marinov et al., 2006;Primeau et al., 2013). However, these modeling studies needed to assume an approximate location of the SOBD to evaluate its relevance. A recent physical oceanography study has identified the surface boundary between the upper and lower overturning cells based on the MOC. Pellichero et al. (2018) explored the MOC in the Southern Ocean in a water mass transformation framework using observations and identified the Neutral density γ = 27.6 kg m −3 outcrop as the dividing line between the upper and lower cells.
Our study builds upon this prior work and aims to constrain the location of the SOBD in greater detail using Lagrangian particle tracking. We released virtual particles throughout the surface Southern Ocean in an eddy-permitting ocean sea-ice model and monitor their trajectories for 15 years when the vast majority of particles have entered either the upper or lower cell of the MOC. This novel approach provides two-dimensional maps of the SOBD as our key outcome. We discuss the implications of the study and compare the SOBD location to observed oceanographic features.

Model and Lagrangian Experiment
This research uses the Australian Community Climate and Earth System Simulator (ACCESS)-OM2-01 model. This ocean-sea ice model has nominally 0.1° horizontal resolution with 75 vertical levels (Kiss et al., 2020). ACCESS-OM2-01 has a good representation of water mass transformation and AABW formation processes when compared to observations (Moorman et al., 2020;Morrison et al., 2020). The model is forced by a JRA55-do repeat year neutral state atmospheric forcing that runs from 1990 to 1991 (Stewart et al., 2020). It was spun up for 50 years, and then daily averaged output saved for 20 years, of which we use the first 16 years in our research (1 year for releasing and the remaining 15 years for running; see below). This model does not include subgrid-scale parameterization for mesoscale eddies since its resolution is considered eddy-resolving, and simulates diffusive mixing in the mixed layer with the K-profile parameterization nonlocal mixing scheme (Kiss et al., 2020). ACCESS-OM2-01 uses a multidimensional piecewise parabolic scheme for advection, with a monotonicity-preserving flux limiter (Kiss et al., 2020).
For the particle tracking experiment we used the Connectivity Modeling System (CMS) v2.0 platform (Paris et al., 2017). CMS uses a fourth-order Runge-Kutta scheme (Paris et al., 2017) to simulate particle trajectories using daily velocity output from the ACCESS-OM2-01 model in the domain south of 40°S. Diffusive mixing processes that impact the movement of tracers in the model are not explicitly included in the model velocity fields used for the particle advection; therefore these diffusive processes are not represented in our particle tracking experiment. Virtual particles were released with a zonally uniform but meridionally enhanced density (see below) 1° × 1° horizontal grid spacing throughout the entire domain (40°S to the Antarctic continent) at 5 and 500 m depth. This release was repeated in four seasons (i.e., the first day of March, June, September, and December) for the first year of the experiment, and then advected for 15 years following the release year (but note that the stream function (Text S1 in Supporting Information S1) during the simulation in the latitude-depth profile over the 16 years. (c) Mean MOC stream function (Text S1 in Supporting Information S1) during simulation in the latitude-2,000 m referenced potential density (σ 2 ) profile over the 16 years. Crosses identify the maximum flux cores of the lower and upper cell and arrows indicate pathways of major water masses. The green-dashed line shows the zonal and time averaged density of the sea surface over the 16 years. Panel (b and c) share the same color bar in (c). In the three plots, red and blue denote the upper and lower cells, respectively. season of the release had no influence on the outcome of the study, which is why we will not further discuss this aspect). Virtual particles' trajectories had a temporal resolution of 1.5 hr and their positions were saved every 5 days, together with the temperature and salinity data at the given locations. Particles reaching topography were reflected and continued the trajectory without interference. The total number of particles released at 5 m (i.e.,"surface") and 500 m was 82,736 and 76,320, respectively (note that the release number was lower at 500 m due to topography). The total number of particles released per 1° × 1° cell during the release year varies from 4 near 40°S to 12 near the Antarctic coast, which corresponds to the finer meridional grid resolution in south since the Mercator grid within ACCESS-OM2-01 (Kiss et al., 2020).
We used density criteria to classify if particles released at any given location in the surface, or 500 m, enter the upper or lower cells. This spatial classification was used to generate maps where the percentage represents the probabilities of particles at any given location south of 40°S to enter the upper or lower cells, to be "not in any cell" or "in both cells". The density criteria (in 2,000 m referenced potential density, σ 2 ) were determined based on the MOC stream function. The MOC stream function in the latitude versus the potential density field (with contours in Sverdrups) shows the existence of two major overturning cells in the Southern Ocean: one is centered at 36.19 kg m −3 , 47°S in latitude, σ 2 coordinates ( Figure 1c, red x), which rotates clockwise and corresponds to the upper cell; another is centered at 37.06 kg m −3 , 65°S, which rotates anticlockwise and corresponds to the lower cell ( Figure 1c, blue x). These two fixed σ 2 density thresholds were used as criteria for the upper and lower cells but have been transformed to surface referenced potential density σ 0 to accurately represent the particle density in the upper ocean. After transformation, the upper cell threshold is σ 0 = 27.15 kg m −3, and the lower cell threshold is σ 0 = 27.80 kg m −3 (detailed information of transformation is provided in the Text S1 and Figure S1 in Supporting Information S1). We note that the nonlinearity of the equation of state leads to nonphysical inversions in potential density in localized regions close to Antarctica, and may have an impact on the location of crossing the lower cell density threshold.
Particles that remain in seawater with σ 0 less than 27.15 kg m −3 for more than 2 years consecutive duration are classified as the upper cell. Upper cell particles account for 77% and 36% in the surface and 500 m releases, respectively. Particles that remain >2 years denser than σ 0 = 27.80 kg m −3 are classified as lower cell and account for 16% and 33% in the surface and 500 m releases, respectively. The 2 years duration was chosen after sensitivity tests have shown that this is a robust criterion for particles to belong to one of the two cells (as discussed in Section 5). Some particles moved north of the northern boundary of the experiment (40°S) before they could satisfy the upper or lower cell criteria (i.e., before crossing for 2 years the upper or lower cell density thresholds). Particles which did this crossing between 0 and 2,000 m depth were classified as upper cell particles as 0-2,000 m falls within the depth range of SAMW and AAIW (Figures 2a and 2c) while particles crossing the boundary below 2,000 m were not classified (i.e., they were counted as "not in any cell" particles). Likewise, particles which remained in the domain for the entire 15 years of the simulation but crossed neither of the density thresholds for >2 years duration were classified as "not in any cell" particles (Figures 2b and 2d). Particles that crossed both thresholds for >2 years duration at some time during the simulation were classified as "in both cells".

SOBD Location
The probabilities of particles entering the upper or lower cell were calculated by counting how many particles released in a 1° × 1° longitude bin eventually satisfy the upper or lower cell criteria. We normalized these counts to the total number of particles released within a bin to get a percentage between 0% and 100%. To present these percentages in one map we subtracted upper cell probabilities from the lower cell probabilities in each bin ( Figure 3). Thus, our definition of the SOBD is the 50% upper cell and 50% lower cell light gray space between the upper (red) and lower cell (blue) in Figure 3. In the following we will use the terms SOBD surf and SOBD 500 , when referring to the 50/50 threshold in the surface ocean or at 500 m depth, respectively.
The SOBD surf is a relatively sharp boundary with a narrow transition from the upper to the lower cell close to the Antarctic continental shelf in most regions (Figure 3a). The sharp boundary is established early in the simulation (after 5-10 years), and there are only 5% of particles that have not crossed the upper or lower cell thresholds after 15 years ( Figure S2i in Supporting Information S1). This indicates that the 15 years of simulation we used in the study is sufficiently long to localize the SOBD surf .
In contrast, the SOBD 500 is a wider circumpolar band with relatively large gray areas ( Figure 3b) where particles (a) have not yet crossed the upper or lower cell thresholds (they are "not in any cell"), or (b) have already crossed both thresholds (they occur "in both cells"). Alternatively (c), there is an equal split between the upper and lower cell. Having a closer look we find strong evidence for the first possibility, since a large fraction from within this wide circumpolar band are "not in any cell" at the end of the 15 years simulation ( Figure S3i in Supporting Information S1). This means that our 15 years simulation is likely not long enough to further narrow down the SOBD 500 . However, we argue that most of these "not in any cell" particles would eventually become part of the upper cell. This is supported by the temporal evolution of the "not in any cell" particles within this circumpolar band, which clearly transition from being "not in any cell" to the upper cell during the 15 years simulation ( Figures S3a-S3c and S3g-S3i in Supporting Information S1). Arguably, most of the gray and slightly red areas in Figure 3b would have become fully red areas if we had let the simulation run longer than 15 years.
Our method of assigning the upper and lower cells based on the defined threshold criteria also leads to some particles being part of both of the cells at different times over the course of the simulation (Figures S2 and S3j-S3l in Supporting Information S1). These particles were not considered for the localization of the SOBD in Figure 3 because they account only for 0.15% and 0.3% of the surface and 500 m depth released particles and it is unclear what their ultimate classification would have been if we had run the simulation for >15 years.

Oceanographic Features Determining the Location of the SOBD
After localizing the SOBD in the model, we checked if the location is correlated with readily observable oceanographic features. This is useful because being able to understand the factors that influence the SOBD location helps to transfer our findings to other models and real-world observations. We specifically looked at surface density outcrops, ACC fronts, subpolar gyre circulation, Ekman transport, and surface buoyancy flux.   and d) show the SOBD 500 . The colorscale represents the percentages of particles in the upper cell minus the percentage of particles in the lower cell at each particle release location. In (a and b), the black line marks the zero meridional Ekman transport contour, the turquoise contour marks the neutral density outcrop γ = 27.6 kg m −3 , which was calculated from the mean density within mixed layer depth as described by Pellichero et al. (2018) using the latest data in Sallée et al. (2021), and the white contours indicate Antarctic Circumpolar Current fronts. From north to south: the Subantarctic Front (Sokolov & Rintoul, 2009); the Southern Boundary (Orsi et al., 1995). In (c and d), the blue and red contours show the σ 0 = 27.15 and 27.80 kg m −3 contours at the surface and 500 m depth. In (b and d), topography shallower than 500 m depth is shaded in light gray. Figure 3 (a and b) show the zero time-mean meridional Ekman transport contour as a black line. The contour was calculated from the wind stress field, which is based on the ACCESS-OM2-01 repeat year wind forcing field (JRA-55-do). The Ekman transport is northward north of the zero meridional Ekman transport contour and southward south of it. The Ekman transport seems to be strongly associated with the location of the SOBD surf throughout Eastern Antarctica and most of the Antarctic Peninsula but shows poorer correlation with it in the Weddell and Ross gyres (Figure 3a). The meridional Ekman transport is weak in a broader region north and south of the zero contour in the Ross and Weddell Seas (not shown), and thus the zero contour is not well defined in these locations. This may explain part of the poorer agreement between the SOBD and zero Ekman transport in these locations. The correlation between the SOBD and zero Ekman transport contour supports earlier modeling studies investigating the biogeochemical relevance of the SOBD, which arbitrarily (or intuitively) defined the SOBD surf location based on Ekman transport Holzer et al., 2014;Primeau et al., 2013). However, we note that the Ekman transport contour is almost identical to the neutral density γ = 27.6 kg m −3 outcrop (which has been identified to demarcate the upper and lower cell by Pellichero et al. (2018)) and the Southern Boundary Front of the ACC in Eastern Antarctica between 30°E and 150°E (Figure 3a). Hence, the SOBD could also be influenced by these two factors within this longitudinal range.
The deviation of the SOBD surf from the Ekman transport contour in the Weddell and Ross gyres may also be influenced by the subpolar gyre circulation. Noticeably, the SOBD surf extends northward to the Southern Boundary in the Ross Sea gyre at around 150°W and to some extent (but much less pronounced) in the Weddell gyre roughly between 0° and 30°E. The gyre circulation enables southward movement of particles released in the surface on the eastern side of the gyres despite the net northward Ekman transport. Thus, the gyre circulation increases the chance for particles released relatively far north to be transported toward the coast and become entrained in forming bottom waters of the lower cell. However, it is important to emphasize that most of particles released in these identified gyre regions still have a higher chance of entering the upper cell than the lower cell (i.e., these regions in Figure 3a are more red than blue).
In addition to the wind-driven circulation, surface buoyancy fluxes ( Figure S4 in Supporting Information S1) influence the location where particles cross the upper and lower cell density thresholds, thus influencing the position of the SOBD. Figure S4 in Supporting Information S1 shows the ACCESS-OM2-01 16-year mean buoyancy heat-equivalent flux. There is not a clear correspondence between the mean buoyancy flux and the SOBD, largely because the particles' density change are influenced by instantaneous fluxes. For example, most of the crossing of particles into the lower cell occurs in winter during AABW formation on the Antarctic continental shelf ( Figure S5 in Supporting Information S1). Differences in buoyancy gain into the upper cell in the Ross and Weddell Seas ( Figure S4 in Supporting Information S1), particularly due to sea ice melt (Abernathey et al., 2016;Pellichero et al., 2018), may contribute to the differences between the Ekman transport and the SOBD in the Ross and Weddell Seas.
The SOBD 500 extends substantially further north than the SOBD surf (Figure 3). This can be explained by the tilted isopycnals (in the depth vs. latitude space) along which the particles upwell (Marshall & Speer, 2012). The most striking result is the close correlation of the edge of 100% lower cell area (i.e., the southern estimate of SOBD 500 ) with the γ = 27.6 kg m −3 outcrop, which has previously been identified by Pellichero et al. (2018) as the boundary between the upper and lower cells in the surface (and thus implicitly as SOBD surf ). Therefore, the good fit of the γ = 27.6 kg m −3 outcrop with the SOBD 500 instead of the SOBD surf as identified in our study was somewhat surprising. This discrepancy may be because a release of virtual particles at 500 m is still shallow enough to allow mixed layer and surface processes to influence particle trajectories. However, it is too deep to be influenced by surface Ekman transport, and the analysis by Pellichero et al. (2018) is based on buoyancy fluxes and not directly on wind patterns.
It is important to note that the γ = 27.6 kg m −3 outcrop largely overlaps with the Southern Boundary of the ACC as determined by Orsi et al. (1995) (Figure 3b). Both of these features represent quite good boundaries of the area belonging 100% to the lower cell (blue area in Figure 3b). The Subantarctic Front corresponds well with the boundary where particles belong 100% to the upper cell (red area in Figure 3b). Thus, the SOBD 500 lies within the ACC boundaries. As mentioned in Section 3 we predict that the SOBD 500 would be narrower and closer to the ACC Southern Boundary than the Subantarctic Front, if the simulation had been extended beyond 15 years. Indeed, the entrainment of some particles into the upper cell may be slow within the ACC because in the ocean interior particles are not exposed to strong surface buoyancy forcing and can therefore move isopycnally for an extended period of time.

Discussion and Conclusions
We defined the SOBD by determining if particles released at the surface, or 500 m depth entered the upper or lower cells based on fixed density thresholds. This raises the question of how sensitive our SOBD maps are to the choice of the density threshold. Perhaps the most convincing argument that our classification is correctly separating the upper and lower overturning cells comes from the individual trajectories of the particles. At both release depths, trajectories of particles classified to the upper or lower cell follow the typical pathway of the Southern Ocean water masses (Figure 2). Upper cell particles moved northward in shallow water, subducting below the mixed layer within SAMW and AAIW. Lower cell particle trajectories generally sink to the abyss and move northward following the path of AABW (Figure 2, Movie S1). Furthermore, the locations where the majority of particles cross the lower cell threshold closely matches deepwater formation regions determined by Morrison et al. (2020) (Figure S5 in Supporting Information S1). Additional sensitivity tests where we altered the threshold densities to minimum/maximum plausible values (i.e., within σ 0 = 27.5 kg m −3 / σ 0 = 27.7 kg m −3 according to the MOC stream function mapped to σ 0 , Figures S1c and S1d in Supporting Information S1) further confirmed that the SOBD pattern is not sensitive to small changes in the density threshold ( Figure S6 in Supporting Information S1). Altogether, these lines of evidence support the notion that our classification scheme is robust.
Another uncertainty in the classification comes from the criterion that a particle must cross a density threshold for at least 2 years to be classified as part of the upper or lower cell. The purpose of this criterion is to reduce noise introduced by short-term diapycnal movement. Indeed, removing the minimum duration limit leads to relatively high percentages of "in both cells" particles (18.9% for the surface and 7.1% for the 500 m release). Already an increase of the limit from 0 to 100 days greatly reduces this problem ("in both cells" particles decreased to 7.8% and 1.9% for the surface and 500 m releases, respectively; Figure S6 in Supporting Information S1). However, to make results more robust, and avoid miscategorizing particles, we chose a relatively long 2 years duration ( Figure  S7 in Supporting Information S1). We also tested the SOBD distribution when only the density of particles at the very end of the 15 years simulation are considered and found that it is similar to the SOBD based on the 2 years minimum duration ( Figure S8c in Supporting Information S1). However, compared with solely using the particles' final density ( Figure S8c in Supporting Information S1), the result using the entire trajectory (Figures 3 and S8a in Supporting Information S1) is more robust since it is less dependent on high frequency diapycnal movement at the end of the simulation. To further test the dependence of the SOBD location on the duration criterion we did additional sensitivity tests with stricter or looser criteria. None changed the appearance of the SOBD noticeably (Figures S8b, S8d and S8e in Supporting Information S1), indicating that our results are robust. As in every model-based assessment, there is uncertainty of the outcome due to the limitations of the ocean model and, in our case, also due to the implementation of virtual particle tracking. Although the ACCESS-OM2-01 model produces dense water formation in the Southern Ocean very well (Moorman et al., 2020), limited model spatial resolution and lack of ice shelf cavities means that there are likely errors in the model representation of the circulation on the Antarctic continental shelf. Furthermore, the lack of explicit mixed layer convection and other interior vertical mixing processes in the particle release experiment is a significant limitation, as without diffusive mixing, the particle trajectories do not fully represent the transport a tracer (i.e., nutrients) would take in the model. Running online model tracer releases would be valuable but requires much larger computational resources. We predict that with the addition of diffusive mixing processes particles would likely cross the upper or lower cell density thresholds more rapidly, and the pattern of the SOBD may differ slightly.
We simulated 15 years with a repeat-year-forcing, that is, the atmospheric forcing of the year 1990-1991 was repeated 15 times. Therefore, our study lacks interannual variability and the influence this may have on the SOBD location. This choice was made (a) because 1990-1991 was an "average year" with no strong influence of major climatic events such as El-Niño Southern Oscillation and Southern Annular Mode (Stewart et al., 2020) and (b) it eliminates the problem to choose in which year to initiate the particle release. Future studies should investigate the sensitivity of the SOBD location to such interannual variability and future climate forcing.
The scientific value of localizing the SOBD can be understood when considering that photosynthetic nutrient consumption north of the SOBD strongly affects surface ocean biogeochemistry outside the Southern Ocean, whereas little such impact occurs when nutrients are consumed south of the SOBD. Several biogeochemical modeling studies have assessed the outstanding significance such spatial differences in nutrient consumption have on global patterns of marine primary production, atmospheric CO 2 draw down, and even plankton community composition (Hauck et al., 2018;Marinov et al., 2006;Nissen et al., 2021;Primeau et al., 2013;Sarmiento et al., 2004). Our study makes the SOBD concept more accessible for the interpretation of real-world data by narrowing down its location. As such, we provide important knowledge that is needed to predict how changes in the spatial distribution of nutrient consumption in the Southern Ocean could affect northern ocean regions. Although uncertainties remain, this constitutes a step forward in our ability to link large-scale physical processes with global biogeochemical element cycling.

Data Availability Statement
The Lagrangian particle trajectory output used in this analysis and the derived data of the SOBD location (as presented in Figure 3) can be found on Zenodo under doi: 10.5281/zenodo.5967481.