The Driving Mechanisms on Southern Ocean Upwelling Change during the Last Deglaciation

: We explore the change in Southern Ocean upwelling during the last deglaciation, based on proxy records and a transient climate model simulation. Our analyses suggest that, beyond a conventional mechanism of the Southern Hemisphere westerlies shift, Southern Ocean upwelling is strongly inﬂuenced by surface buoyancy forcing and the local topography. Over the Antarctic Circumpolar Current region, the zonal mean and local upwelled ﬂows exhibited distinct evolution patterns during the last deglaciation, since they are driven by different mechanisms. The zonal mean upwelling is primarily driven by surface wind stress via zonal mean Ekman pumping, whereas local upwelling is driven by both wind and buoyancy forcing, and is tightly coupled to local topography. During the early stage of the last deglaciation, the vertical extension of the upwelled ﬂows increased downstream of submarine ridges but decreased upstream, which led to enhanced and diminished local upwelling, downstream and upstream of the submarine ridges, respectively.


Introduction
Paleo-proxy reconstructions have shown dramatic changes in the Southern Ocean upwelling during the last deglaciation, which potentially contribute to the deglacial rise of atmospheric CO 2 [1][2][3][4][5]. The "westerlies shift" hypothesis is one of the mostly discussed mechanisms on the changes in Southern Ocean upwelling. According to this hypothesis, an equatorward shift of the Southern Hemisphere westerlies acts to weaken the Southern Ocean upwelling during glacial periods such that more respired CO 2 accumulates in deep oceans while the atmospheric CO 2 concentration reduces [6][7][8]. During interglacial periods, a poleward displacement of the the Southern Hemisphere westerlies intensifies the Southern Ocean upwelling and respired CO 2 outgassing [1,9,10], thus leading to an increase in CO 2 in the atmosphere. Such glacial-interglacial swings of the Southern Hemisphere westerlies and Southern Ocean upwelling could be coupled to atmospheric CO 2 variations and be related to the global mean climate, the bipolar seesaw and Earth's obliquity [5].
The "westerlies shift" hypothesis is intrinsically proposed for the upwelled flow over the Southern Ocean from a zonal mean perspective, whereas evidences of Southern Ocean upwelling change are mostly derived from a limited number of proxy records at scattered sites [2,31,32]. As shown in previous studies, perhaps due to the lack of an integrated data set, the local upwelling from individual proxy record is usually treated as a representative of zonal mean upwelling, and therefore, the wind-driven mechanism is adopted to explain the glacial-interglacial changes in Southern Ocean upwelling from a zonal mean perspective [2]. However, given the complex local topography and, in turn, the significant differences between local and zonal mean upwelled flows, it has remained unclear whether features derived from a few local upwelling proxies can truly represent those of zonal mean upwelling over the Southern Ocean.
In this study, we elucidate that local and zonal mean upwelled flows in the Southern Ocean are driven by distinct mechanisms. Combining paleo-proxy records and a simulation of the Transient Climate of the Last 21,000 Years (TraCE-21ka) [33], we show that, besides the "westerlies shift" hypothesis, several alternative mechanisms, such as freshwater forcing and topography effects, can explain the changes in Southern Ocean upwelling during the last deglaciation. Our mechanisms are built on Southern Ocean dynamics [34,35] in which the bottom topography in the Southern Ocean can strongly affect the wind-driven upwelling by undermining the barotropic structure of vertical flows. This topographic effect, meanwhile, is largely dependent on deep ocean stratification. During the last deglaciation, deep ocean stratification was greatly modulated by the changes in surface freshwater fluxes over the Southern Ocean. As a result, topographic effects expectedly varied correspondingly during the deglaciation, and should have contributed to the change in Southern Ocean upwelling.
The rest of the paper is structured as follows: Section 2 gives a brief introduction of TraCE-21ka and paleo-proxy data as well as the approaches used for the data-model comparison. A subsequent analysis is made in Section 3 that focuses on the mechanisms driving the changes in Southern Ocean upwelling during the last deglaciation from local and zonal mean perspectives. Conclusions and discussions are given in Section 4.

Model and Approaches
We use a synchronously coupled atmosphere-ocean general circulation model, the National Center for Atmospheric Research (NCAR) Community Climate System Model version 3 (CCSM3) [36], in the version of T31_gx3v5 resolution [37]. CCSM3 includes the atmosphere component of Community Atmosphere Model version 3 (CAM3) [38], the land component of Community Land Surface Model version 3 (CLM3) [39], the sea ice component of Community Sea Ice Model version 5 (CSIM5) [40] and the ocean component of Parallel Ocean Program version (POP) [41]. CAM3 and CLM3 adopt a horizontal resolution of T31 spectral truncation (3. TraCE-21ka is carried out with CCSM3 and initialized from an earlier LGM equilibrium simulation in which dynamic vegetation is included [42]. From LGM to the present-day climate, TraCE-21ka adopts the boundary conditions and external forcings, such as Earth's orbits, atmospheric greenhouse gases (GHGs), ice sheets and meltwater discharge into oceans inferred from paleo-evidences (Figure 1a,b). Driven by these boundary conditions and external forcings, TraCE-21ka is successful in capturing many key features of climate evolution during the last deglaciation, such as changes in the Atlantic Meridional Overturning Circulation (AMOC; Figure 1c) as well as cooling during Heinrich event 1, the abrupt Bølling-Allerød warming and cooling in the Younger Dryas [33,43,44]. In this study, analyses are performed using decadal mean model outputs of TraCE-21ka during the last 22,000 years.  [45] and atmospheric CO 2 concentration (green) [46]; ppmv, parts per million by volume. (b) Freshwater fluxes (FWF) in the Northern (dark pink) and Southern (light pink) Hemispheres in TraCE-21ka. (c) The 231 Pa/ 230 Th ratio (brown) at Bermuda (GGC5 core) as a proxy for AMOC export [47], and model maximum AMOC transport (below 500 m, orange). (d-h) The column integrated upwelling (gray) and its 100-year average (black) in TraCE-21ka at corresponding grids of proxy cores (d) A, (e) B, (f) C and (g) D, and (h) averaged over the ACC region. In these panels, opal fluxes (blue) and initial unsupported 231 Pa/ 230 Th ratio (red) are reproduced from previous studies [2,32]. In (h), the maximum of the 100-year mean Southern Hemisphere westerlies on 867 hpa is shown in orchid. As the shift of the westerlies could be small and within a grid of the atmosphere model, a quadratic interpolation around the grid point with maximum mean wind is adopted to yield the exact latitude of the mean jet [25]. H1: Heinrich event 1; BA: Bølling-Allerød; YD: Younger Dryas.
Regarding observations, we use preserved opal fluxes and the 231 Pa/ 230 Th ratio in the sedimentary record as proxies for Southern Ocean upwelling [2,32]. The benchmark here is to compare the simulated upwelling with its proxy-indicated counterpart at multiple sites for the credibility of further analyses. Since TraCE-21ka does not include a silica cycle in its ocean component, we cannot make a direct comparison between the model and observations in terms of opal fluxes or other upwelling proxies. As an alternative, we use the column integrated upwelling as a proxy analogous to opal fluxes, considering that opal production is likely to be governed by an integrated upwelled flow throughout different ocean depths [48]. At individual sites, the local upwelling is estimated as a gridded upward velocity w up on each level and a gridded volume-averaged velocity w up over the whole depth. Particularly, w up and w up are calculated as follows: where w up (k) and dz(k) are, respectively, the gridded upwelling velocity and the thickness of the kth level, and N is the number of grids in the vertical at each site. Moreover, to explore the Southern Ocean upwelling from a zonal mean perspective, we focus on zonal mean vertical velocity w and Ekman pumping w Ek . We examine the area-averaged velocity (w ACC ) over the latitudinal belt of the Antarctic Circumpolar Current (ACC) and the upwelling over the ACC region (w up ACC ). Particularly, w ACC and w up ACC are calculated as follows: where w up ACC (k) is the value of w up ACC at the kth level. y N and y S are the northern and southern boundaries of the ACC region and A xy denotes the area of the ACC region. We also calculate the Ekman pumping averaged over the ACC region: denotes the difference of zonal wind stress between the southern (φ S = 62 o S) and northern (φ S = 55 o S) boundaries, w e denotes local Ekman pumping, ∆φ = φ N − φ S denotes the latitudinal span, f is Coriolis parameter and ρ 0 is water density.

Southern Ocean Upwelling Indicated by Proxy Records and in TraCE-21ka
We investigate the local upwelling indicated by paleo-proxies and simulated by TraCE-21ka at four sites ( . These cores are selected to provide high-resolution records of opal accumulation during the last deglaciation, which are denoted as cores A-D, thereafter. Since the LGM, upwelling implied from proxy records exhibited dramatic changes. It increased toward ∼15 ka at cores A, B and C, but decreased at core D ( Figure  1d-g), which indicates that changes in zonal mean and local upwelled flows over the ACC region potentially followed different evolution patterns during the last deglaciation. The simulated upwelling in TraCE-21ka, on the other hand, exhibits generally consistent features with the proxy records over the four sites, albeit with some differences in the timing and magnitude of the upwelling changes (Figure 1d-g). Particularly, both observations and model simulations reveal intensified upwelling at cores A-C but weakened upwelling at core D. In contrast, the zonal mean upwelling over the ACC region has little altered during the last 22,000 years, except a small increase during the Bølling-Allerød interstadial (Figure 1h), which is different from that of any local upwelling at cores A-D. Moreover, it is worth noting that TraCE-21ka simulates a poleward displacement of Southern Hemisphere westerlies during the LGM from their present day position (Figure 1h). This pattern is consistent with the wind shifts in multiple PMIP3 models, indicative of the role of Antarctic sea ice in modulating the westerlies shift [28]. During the last deglaciation, Southern Hemisphere westerlies showed a consistently equatorward shift in TraCE-21ka (Figure 1h), suggesting that the westerlies shift seems unlikely to be the major mechanism for the change in Southern Ocean upwelling.
Based on the results above, we elucidate that the available local upwelling reconstructions from individual proxy records are not representative of the zonal mean upwelling over the ACC region and thus cannot be interpreted to be forced predominantly by the westerlies shift. Instead, changes in the local and zonal mean ocean upwelled flows during the last deglaciation could be driven by different mechanisms.

The Mechanism of Zonal Mean Upwelling Changes in the Southern Ocean
We find that the change in zonal mean Southern Ocean upwelling is driven mainly by changes in surface wind stress and associated Ekman pumping. As shown in Figure 3a, surface westerly winds over the Southern Ocean drive a northward Ekman transport, whose divergence leads to a positive (upward) Ekman pumping to the south of the wind maximum, and in turn, an upwelling below the Ekman layer (Figure 3c). Over the ACC region, the wind effect can penetrate to a mid-depth of ∼2000 m, above which the ocean is unblocked by the topography and is, therefore, driven by interfacial form stress, according to the "open channel" dynamics [49,50]. As a result, the magnitude of the zonal mean upwelling in the upper 2000 m equals approximately that of the zonally mean Ekman pumping within this region (Figure 3d). Note here that the definition of "approximate equality" is based on two conditions: (1) the correlation coefficient between vertical velocity w and Ekman pumping w Ek is higher than 0.9 (see Figure 3b); and (2) the difference between the two is one order smaller than either in magnitude (see Figure 3d). During the last deglaciation, the zonal mean Ekman pumping and therefore, the vertical velocity in the upper ∼ 2000 m over the ACC showed a general decline trend, except a sharp increase during the Bølling-Allerød interstadial (Figure 4e). Below ∼2000 m, waters are blocked by the Drake Passage and flows are greatly affected by the bottom topography. Surface wind forcing can penetrate the unblocked levels over the ACC region and drive deep ocean flows in an advective ventilation depth [51]. As a result, the deep ocean upwelling extends deepest during the Bølling-Allerød interstadial (Figure 4e). To summarize, our results suggest that the zonal mean upwelling is primarily wind-driven in the ACC region, exhibiting a decline trend throughout the last 22,000 years, except a short intensification during the Bølling-Allerød interstadial.

The Mechanism of Local Ocean Upwelling Changes in the Southern Ocean
We further delve into the mechanism that determines the local upwelling at the four core sites, using TraCE-21ka. We find that changes in the column integrated upwelling (Figure 1d-g) depend more on the vertical extension than the strength of upwelled flows during the last deglaciation. This feature is evident in the vertical profiles of local upwelling at the four core sites (Figure 4a-d). From ∼17 ka to ∼14 ka, increases of local upwelling at sites A, B and C are mainly caused by a vertical extension of upwelled flows from ∼1000 m to ocean bottom (Figure 4a-c), whilst the decrease of the local upwelling at site D is primarily related to shallowing upwelled flows from the ocean bottom to ∼1000 m ( Figure  4d). These changes in the vertical extension of local upwelled flows are much larger than that in zonal mean upwelled flows over the ACC region (Figure 4e).
The difference in the change of vertical extension between local and zonal mean upwelled flows indicates different mechanisms driving these flows during the last deglaciation. Herein we examine the correlation between vertical velocity and Ekman pumping from both local and zonal mean perspectives. We find that the zonal mean vertical velocity over the ACC region is highly correlated with the zonal mean Ekman pumping to a depth of ∼2000 m in the unblocked levels ( Figure 5d). In contrast, correlation profiles between local Ekman pumping and local vertical velocity vary greatly along a latitudinal circle (Figure 5a,b). The zonal average of these local correlations diminishes rapidly with depth ( Figure 5c,d), indicating that local upwelling is driven by mechanisms beyond local wind forcing. To better understand the mechanism that controls the vertical extension of local upwelling, we interpret the vertical velocity w into two parts: the local Ekman pumping w e , and the depth-dependent w pva that denotes the vertical integration of the convergence of the meridional velocity of large-scale geostrophic flows. Particularly, if we start from the planetary vorticity conservation in a water column below the Ekman layer, we have the following: where β is the gradient of the Coriolis parameter f , and v denotes the meridional velocity.
Integrating Equation (7) from an arbitrary depth z to the bottom of the Ekman layer z = −d, we have the following: where the Ekman pumping w e = w(−d) = curl( τ 0 ) , τ 0 is the surface wind stress, and where the integration is from an arbitrary depth z that is away from ocean bottom to the bottom of the Ekman layer. Here, w pva is related to a local topographic steering [45]. Based on Equation (8), the vertical velocity, and in turn, the local upwelling, are expected to be different downstream (upstream) of a submarine ridge where a strong southward (northward) flow exists. During the LGM, cores A, B and C were located downstream of the Macquarie Ridge, the Kerguelen Plateau, and the American-Antarctic Ridge, respectively (Figure 2), where a prevailing southward flow (Figure 6a-c) was accompanied by a vertical stretching that increased with depth. As a result, the upward velocity was reduced from the base of the Ekman layer with depth, and may have vanished or even been reversed to downwelling at certain depths. In contrast, core D was located in the upstream of the Macquarie Ridge (Figure 2), so the prevailing northward flow there (Figure 6d) was accompanied by a compression of the water column, which enhanced the upwelling from the Ekman layer toward deeper levels. Clearly, meridional vorticity advection controls the vertical extension of upwelled flows below the Ekman layer and therefore, regulates the variation of local upwelling. Figure 4f-i shows the decomposition of the vertical velocity w into w e and w pva at core sites A-D during the past 22,000 years. It is evident that changes in the vertical extension of local upwelling at four sites are determined primarily by w pva . Furthermore, the change of w pva is associated with the vertical shear of the meridional flow and therefore, the zonal density gradient, according to the thermal wind relationship. Since the meridional flow is approximately geostrophic in places away from the Ekman layer and ocean bottom, the thermal wind relationship can be written as follows: ∂v ∂z where v g denotes the geostrophic meridional velocity, g is the gravitational acceleration and ρ denotes the local water density. Figure 7 displays the thermal wind relationship at the four cores, which relates a negative (positive) shear of the meridional flow at cores A-C (core D) to the eastward decreasing (increasing) density in the Southern Ocean. Accordingly, changes in zonal density gradients alter the vertical shears of meridional flows. During the LGM, cores A-C (core D) showed strong negative (positive) zonal density gradients. However, these negative (positive) density gradients reduced in magnitude during the last deglaciation. As a result, the vertical shears of meridional flows at all four cores became smaller in magnitude, correspondingly.
To unravel the reason on the change in zonal density gradients, we write the zonal density gradient (∂ρ/∂x) as follows: where α = − 1 ρ ∂ρ ∂t p,s is the thermal expansion coefficient and β = 1 ρ ∂ρ ∂s p,t is the saline contraction coefficient, and −α ∂t ∂x and β ∂s ∂x denote the thermal and haline contributions to the zonal density gradient, respectively. We find that changes in zonal density gradient are mainly caused by the haline contribution, i.e., due to the changes in zonal salinity gradient at all four sites (Figure 7). These changes in zonal salinity gradient, moreover, are closely tied to the changes in the meridional salinity gradient. Herein we examine the salinity budget as follows: − us x − vs y − ws z + Res = 0 (11) which states a balance among the zonal (−us x ), meridional (−vs y ) and vertical (−ws z ) salinity advection and a residual term (Res) that includes salinity tendency and parameterized sub-grid mixing. We find that, at cores A-C (core D), salinity budgets are dominated by southward, negative (northward, positive) meridional salinity advection and eastward, positive (westward, negative) zonal salinity advection (Figure 8). From the time of Heinrich event 1 toward the Bølling-Allerød interstadial, meridional salinity gradients diminish over the ACC region, which leads to reduced meridional and zonal salinity advection and diminished zonal salinity gradients across the four proxy sites. Furthermore, we find that the reduction in meridional salinity gradient is generated by changes in surface freshwater fluxes over the Southern Ocean. The total surface freshwater flux (SFWF) can be written as follows: where SFWF is composed of various freshwater fluxes from precipitation (P), evaporation (E), river runoff (R), sea ice melt (M) and brine rejection (Br) due to the melting of sea ice. We find that, during the last deglaciation, changes in surface freshwater fluxes over the Southern Ocean were essentially contributed by those in sea ice melt and brine rejection ( Figure 9). From the LGM to early Holocene, surface ocean temperature around Antarctica significantly increased due to a bipolar see-saw response [52] to meltwater discharge into the North Atlantic [53] during Heinrich event 1 and the Younger Dryas as well as a CO 2 warming [54]. Accordingly, less Antarctic sea ice was formed in austral winter such that brine rejection was also reduced, which led to an anomalous freshwater input around Antarctica to dilute Antarctic Bottom Water (AABW). During the last deglaciation, the dilution of AABW was even stronger than that of North Atlantic Deep Water (NADW), which is indicated from the proxy records [55,56]. Therefore, the meridional salinity gradient across the Circumpolar Deep Water (CDW), as exampled by the salinity contrast between AABW and the southward outflow of NADW through the South Atlantic sector, was reduced ( Figure 10).

Conclusions and Discussion
In this study, we combined paleo-proxy observations and a fully coupled climate model simulation, TraCE-21ka, to analyze the change in Southern Ocean upwelling during the last deglaciation. We tested the "westerlies shift" hypothesis but found that the displacement of Southern Hemisphere westerlies cannot explain the upwelling change as indicated from paleo-proxies and TraCE-21ka. We then detected the driving mechanisms on the Southern Ocean upwelling change from both zonal mean and local perspectives. We showed that the zonal mean upwelling over the ACC region is primarily driven by surface wind stress via zonal mean Ekman pumping, whilst the local upwelling is driven by both surface wind stress and surface buoyancy forcing and is tightly coupled to local topography. The latter mechanism can be summarized as follows: dwindling Antarctic sea ice during the last deglaciation helped reduce the meridional salinity gradient in the Southern Ocean, leading to diminished negative (positive) zonal salinity gradients downstream (upstream) of submarine ridges. Due to a dominant haline contribution to density change, the magnitude of zonal density gradient decreased, which abates the baroclinicity in meridional flows and hence associated planetary vorticity advection. As a result, the vertical extension of upwelled flows became deeper (shallower) downstream (upstream) of submarine ridges, corresponding to intensified (weakened) column integrated upwelling.
Our results show that Antarctic sea ice and associated surface freshwater fluxes played a vital role in modulating Southern Ocean upwelling during the last deglaciation, which is consistent with the results from previous studies [57][58][59][60][61][62][63][64][65][66][67] that highlight the effects of buoyancy forcing on Southern Ocean overturning and upwelling. While most of these studies compared ocean circulations between LGM and modern climates, they suggested that the expansion of Antarctic ice extent during the LGM led to a contraction of the overturning of NADW but an expansion of the overturning of AABW, hence pushing the Southern Ocean upwelling equatorward and abating respired CO 2 outgassing. Our results, furthermore, demonstrate that, besides the zonally integrated overturning and upwelling, Antarctic sea ice and associated surface freshwater forcing can also significantly influence the local upwelling from LGM to the modern climate.
Eddy effects are not discussed in this study due to a lack of model outputs of eddyinduced velocities in TraCE-21ka. It was suggested that the Eulerian-mean meridional overturning circulation (MOC) can be partially offset by a reversed eddy-induced MOC, resulting in a residual MOC [68,69] in the Southern Ocean. Meanwhile, as we have shown in this study, the column integrated upwelling in TraCE-21ka can still capture the robust change in Southern Ocean upwelling during the last deglaciation as indicated by proxy records. This is likely due to a major contribution from the Eulerian-mean upwelling to the tracer's movements. Nevertheless, it is worth noting that eddies still play a role in modulating Southern Ocean upwelling. We estimate that oceanic eddies may partially modify the change in the Eulerian mean upwelling either from a local or zonal mean perspective, due to their essences to even the baroclinicity in the Southern Ocean.
It is also worth noting that uncertainties exist in the impacts of Southern Ocean upwelling change on the carbon cycle and atmospheric CO 2 variations. Several climate model simulations depict that changes in Southern Ocean upwelling only weakly influence atmospheric CO 2 variations [70,71], whilst many others report that Southern Ocean upwelling changes could well explain the atmospheric CO 2 changes during the last deglaciation [10,72]. Moreover, besides the Southern Ocean upwelling change, there are many other processes and mechanisms, such as biological pump and vertical mixing, potentially affecting the carbon sequestration in deep ocean and respired CO 2 outgassing [4,66,73].