The Roles of Wind and Sea Ice in Driving the Deglacial Change in the Southern Ocean Upwelling: A Modeling Study

: The Southern Ocean (SO) played a fundamental role in the deglacial climate system by exchanging carbon-rich deep ocean water with the surface. The contribution of the SO’s physical mechanisms toward improving our understanding of SO upwelling’s dynamical changes is devel-oping. Here, we investigated the simulated transient SO atmosphere, ocean, and sea ice evolution during the last deglaciation in a fully coupled Earth system model. Our results showed that decreases in SO upwelling followed the weakening of the Southern Hemisphere surface westerlies, wind stress forcing, and Antarctic sea ice coverage from the Last Glacial Maximum to the Heinrich Stadial 1 and the Younger Dryas. Our results support the idea that the SO upwelling is primarily driven by wind stress forcing. However, during the onset of the Holocene, SO upwelling increased while the strength of the wind stress decreased. The Antarctic sea ice change controlled the salt and freshwater ﬂuxes, ocean density, and buoyancy ﬂux, thereby inﬂuencing the SO’s dynamics. Our study highlighted the dynamic linkage of the Southern Hemisphere westerlies, ocean, and sea ice in the SO’s latitudes. Furthermore, it emphasized that zonal wind stress forcing and buoyancy forcing control by sea ice together regulate the change in the SO upwelling.


Introduction
The last deglaciation from 21 to 9 kyr BP (before present) offers a case study for the most recent natural global warming process. This period involved a transient climate evolution from a colder to a warmer world. It was interrupted by Northern Hemisphere millennial events with a simultaneous CO 2 rise that is comparable to anthropogenic emission since the industrial revolution [1]. The Southern Ocean (SO) plays a pivotal role in controlling the variation in atmospheric CO 2 concentrations over glacial-interglacial cycles [2]. It ventilates a significant part of the global ocean and acts as a window for the return of carbon-rich deep ocean water to the surface. The dynamical theory of the SO dynamics is developing. Thus, understanding the physical mechanisms that determined the temporal and spatial change in SO upwelling throughout the Earth's history will improve our comprehension of the response of the SO's dynamics in future climate projections and the associated difference in atmospheric CO 2 concentrations.
During the northern Atlantic millennial events, studies show a correlation between the change in the atmospheric CO 2 concentrations and Antarctic temperatures involving a poleward shift in the Southern Hemisphere mid-latitude westerlies and an increase in the SO overturning [3]. Northern Atlantic cooling during the Heinrich Stadial 1 and the Younger Dryas overlapped with the Southern Hemisphere deglacial warming, with a rise in wind-driven SO upwelling [4][5][6][7][8][9]. Thus, Northern Hemisphere stadial millennial events,  [11]. (b,f) Northern hemisphere meltwater (NHMW) forcing (black) and Southern Hemisphere meltwater (SHMW) forcing (red, dashed) (units are in meters per kiloyear (m kyr -1 )). Simulated upwelling (black; units are in Sverdrup (Sv)) and the maximum zonal wind stress (≈46° S; red, dashed; units are in dyne per square centimeter (dyne cm -2 )) was smoothed using a Savitsky-Golay filter with a third-order polynomial convolution at the Indian Ocean (IO) (c,g) and the Pacific Ocean (PO) sector (d,h) respectively. The simulations were area-averaged over the IO sector (42°S to 58°S and 86°E to 140°E) and the PO sector (42°S to 58°S and 157°W to 86°W).

Analyzing the Performance of the CCSM3 TraCE-21ka Model Regarding Simulating the Present-Day Southern Hemisphere Climate
This section presents the analysis of the performance of CCSM3 TraCE-21ka when simulating the SO atmospheric and oceanic features in both the Indian Ocean (IO) and Pacific Ocean (PO) sectors with the present-day NCEP/NCAR reanalysis. The model re-  [11].
(b,f) Northern hemisphere meltwater (NHMW) forcing (black) and Southern Hemisphere meltwater (SHMW) forcing (red, dashed) (units are in meters per kiloyear (m kyr −1 )). Simulated upwelling (black; units are in Sverdrup (Sv)) and the maximum zonal wind stress (≈46 • S; red, dashed; units are in dyne per square centimeter (dyne cm −2 )) was smoothed using a Savitsky-Golay filter with a third-order polynomial convolution at the Indian Ocean (IO) (c,g) and the Pacific Ocean (PO) sector (d,h) respectively. The simulations were area-averaged over the IO sector ( The employed coupled atmospheric model was the Community Atmospheric Model 3 (CAM3) performed at a T31 (≈3.75 degrees) resolution with a vertical resolution of 26 hybrid coordinate levels. The coupled oceanic component (Parallel Ocean Program (POP)) is in spherical polar coordinates with a dipole grid. It has a resolution of gx3v5 (about 3.6 degrees longitudinal resolution with a variable latitudinal resolution; there is a finer resolution to about 0.9 degrees near the equator), with a vertical resolution of 25 levels. The sea ice model employed was the NCAR Community Sea Ice Model (CSIM), which includes a subgrid-scale ice thickness distribution. Its resolution is identical to the POP model. The ice and ocean model outputs were gridded to a T31 resolution for consistency with the atmospheric model. The ocean model caveat was that owing to its coarse resolution, it fails to resolve ocean eddies and instead uses the Gent and McWilliams parameterization [47,48]. However, the latest studies have advocated that coarse resolution models agree with the high-resolution models [15] when simulating a realistic Antarctic Circumpolar Current and the SO overturning response to freshwater and wind stress forcing. The model output is publicly available for download from the National Center for Atmospheric Research Climate Data Gateway. TraCE-21ka model simulations have been shown to replicate many vital features of the global climate evolution during the last deglaciation [46,49], including El Niño [50], Asian-African monsoon [51], and high-latitude temperature [52]. It also demonstrated good agreement regarding the simulated atmospheric regional surface air temperature over Antarctica, the Atlantic Meridional Ocean Circulation transport, and local sea surface temperature in the south Pacific and the SO with paleo data [53]. Additionally, TraCE-21ka simulations reproduced the Southern Hemisphere atmospheric circulation's general characteristics well with present-day conditions [45].

Analyzing the Performance of the CCSM3 TraCE-21ka Model Regarding Simulating the Present-Day Southern Hemisphere Climate
This section presents the analysis of the performance of CCSM3 TraCE-21ka when simulating the SO atmospheric and oceanic features in both the Indian Ocean (IO) and Pacific Ocean (PO) sectors with the present-day NCEP/NCAR reanalysis. The model reproduced near-surface 950 hPa zonal wind speeds ( Figure S1) with a similar seasonal pattern compared to the present-day NCEP/NCAR reanalysis [54]. Near-surface Southern Hemisphere westerly winds were strongest in the ASO (August, September, and October) season compared to the JJA (June, July, and August). In the PO sector, a 33-year (1958 to 1990) monthly mean 950 hPa zonal wind climatology had strong spatial (r = 0.8 and R 2 = 0.97) and ASO seasonal correlations (r = 0.99, R 2 = 0.81, p = 0.003, and RMSE = 0.2) with the NCEP reanalysis I. In the IO sector, the simulated near-surface zonal wind had a strong ASO spatial correlation (r = 0.92, R 2 = 0.94, p < 0.05, and RMSE = 0.2) with the NCEP reanalysis I (one month offset). Additionally, 70-year (1920 to 1990) annual average decadal mean sea surface temperature climatology from the TraCE-21ka POP model ( Figure S2) had a strong correlation in the IO (r = 0.9, R 2 = 0.5, and p = 0.006) and PO (r = 0.8, R 2 = 0.5, and p = 0.05) with the National Oceanic and Atmospheric Administration (NOAA) Extended Reconstructed Sea Surface Temperature (SST) V5 [55]. Simulated sea ice and ocean temperatures throughout the SO showed reasonable agreement with the IO and PO sectors' proxy records [56]. TraCE-21ka underestimated the sea surface temperatures. Nonetheless, our evaluations support the idea that the TraCE-21ka model reproduces the general properties of the SO wind and oceanic changes in the present-day reasonably well.

Study Domain
To investigate the last deglacial change in the SO upwelling, we chose four millennial timescales: (1) the Last Glacial Maximum (LGM; 20 to 19 kyr BP), (2) the Heinrich Stadial 1 (H1; 17 to 16 kyr BP), (3) the Younger Dryas (YD; 12.5 to 11.5 kyr BP), and (4) the onset of the Holocene (O_H; 10 to 9 kyr BP). These four periods illustrate Earth's transient climate evolution with changing Southern Hemisphere westerlies and sea ice history. Realistic simulation of the Antarctic Cold Reversal (≈14.5 to ≈12.5 kyr BP) has been Sustainability 2021, 13, 353 5 of 21 a challenge for the modeling community [56]. TraCE-21ka simulates Bølling-Allerød (≈14.5 kyr BP) warming by varying the meltwater forcing with a sudden termination, followed by the abrupt reopening of the Northern Hemisphere and Southern Hemisphere meltwater flux (Meltwater Pulse 1-A). This massive influx of meltwater at the Ross Sea and Weddell Sea surfaces possibly strengthens the upper meridional overturning cell [19], which leads to the SO upwelling's unrealistic overshooting. Therefore, in this study, we excluded this modeling limitation of the SO upwelling's overestimation during the Antarctic Cold Reversal.
Unlike other modeling studies that focused on the SO average zonal perspective [7], we analyzed the IO and PO at the basin scale. To identify appropriate representative sectors, we selected sections of the IO and PO basins that overlapped with the simulated Southern Hemisphere wintertime maximum near-surface winds, as shown in Figure 2. The deepest mixed layer depths are located in the Subantarctic region equatorward to the SO upwelling zone, where the Subantarctic mode water is formed and transported equatorward [57]. These representative sectors are poleward to the simulated annual mean mixed layer depths in the SO upwelling zone (simulated negative wind stress curl), as illustrated in Figure 2. Additionally, the sectors are north of the sea ice boundary and receive freshwater and salt fluxes that are associated with the Antarctic sea ice melting and freezing. Since stronger surface zonal wind may result in stronger Eulerian mean wind-driven overturning, we examined the ASO seasonal atmospheric climatology. A limitation of our ocean and sea ice model analysis is that we used annual mean data instead of seasonal due to data availability. However, the limited available ASO monthly mean ocean model output captured the ASO season's general characteristics in the annual mean ocean output (r ≈ 0.9, p < 0.05). The simulated SO upwelling volume was calculated as an average vertical velocity in the representative sectors integrated to a 500 m depth. We reasoned that the vertical integration to deeper ocean velocities would better represent the SO upwelling. The conclusions drawn from this study remained the same, irrespective of the study area's size and depth, as long as it overlapped with maximum surface winds. LGM, H1, YD, and O_H represent the Last Glacial Maximum, the Heinrich Stadial 1, the Younger Dryas, and the onset of Holocene, respectively. White represents the easterly zonal wind and a lack of model values.

Zonal Wind
Simulation of the near-surface 950 hPa zonal wind indicated that there were strong (10 to 15 m s -1 ) winds throughout the deglaciation. The representative sectors depicted in black boxes overlapped it and lay poleward to 80 m annual mixed-layer depth contours ( Figure 2). The simulation illustrated stronger near-surface wind and deeper yearly mixed-layer depths in the IO than the PO sector at all timescales. In particular, the model shows that the LGM had the strongest near-surface SO zonal wind (0.5 to 1.5 m s -1 higher) as compared to other time slices during the last deglaciation ( Figure 3). LGM, H1, YD, and O_H represent the Last Glacial Maximum, the Heinrich Stadial 1, the Younger Dryas, and the onset of Holocene, respectively. White represents the easterly zonal wind and a lack of model values.

Zonal Wind
Simulation of the near-surface 950 hPa zonal wind indicated that there were strong (10 to 15 m s −1 ) winds throughout the deglaciation. The representative sectors depicted in black boxes overlapped it and lay poleward to 80 m annual mixed-layer depth contours ( Figure 2). The simulation illustrated stronger near-surface wind and deeper yearly mixedlayer depths in the IO than the PO sector at all timescales. In particular, the model shows that the LGM had the strongest near-surface SO zonal wind (0.5 to 1.5 m s −1 higher) as compared to other time slices during the last deglaciation ( Figure 3).   Figure 1d,h show the deglacial SO upwelling in the IO and the PO sectors, respectively. Spatially, the model simulated more vigorous upwelling in the IO than the PO sector at all timescales, which is consistent with present-day ocean circulation. The simulated upwelling core remained intact during all the time slices at both the ocean sectors, inferring that there was no significant spatial change in the deglacial Antarctic Circumpolar Current circulation. The model simulated the most vigorous upwelling during the LGM (≈7 Sv) in the IO sector, which was 5 to 7% stronger than in the H1, YD, and O_H (Figure 1c,g). The simulation showed that the upwelling core (94°E to 109°E) weakened from the LGM to the YD and then reinvigorated in the O_H. In particular, in the PO sector, the O_H showed the greatest upwelling strength among the four millennial timescales ( Figure 1d). Upwelling in the O_H (≈4 Sv) was 4% stronger than in the LGM, indicating more vigorous upwelling at the end of deglaciation. Table 1 presents an overview of it.

Zonal Wind Stress
The simulated wind stress during all the time slices in the IO sector was ≈36% stronger than the PO sector, in concordance with the ≈37% stronger upwelling volume.   Figure 1d,h show the deglacial SO upwelling in the IO and the PO sectors, respectively. Spatially, the model simulated more vigorous upwelling in the IO than the PO sector at all timescales, which is consistent with present-day ocean circulation. The simulated upwelling core remained intact during all the time slices at both the ocean sectors, inferring that there was no significant spatial change in the deglacial Antarctic Circumpolar Current circulation. The model simulated the most vigorous upwelling during the LGM (≈7 Sv) in the IO sector, which was 5 to 7% stronger than in the H1, YD, and O_H ( Figure 1c,g). The simulation showed that the upwelling core (94 • E to 109 • E) weakened from the LGM to the YD and then reinvigorated in the O_H. In particular, in the PO sector, the O_H showed the greatest upwelling strength among the four millennial timescales ( Figure 1d). Upwelling in the O_H (≈4 Sv) was 4% stronger than in the LGM, indicating more vigorous upwelling at the end of deglaciation. Table 1 presents an overview of it.

Zonal Wind Stress
The simulated wind stress during all the time slices in the IO sector was ≈36% stronger than the PO sector, in concordance with the ≈37% stronger upwelling volume. The simulated upwelling core overlapped spatially with the wind stress maxima. The SO upwelling in the IO sector had a strong correlation with wind stress maxima (46 • S) (r ≈ 0.8; p < 0.05) (Figure 1c,g), whereas, in the PO sector, it had a Pearson correlation coefficient of about 0.2 (p < 0.05) (Figure 1d,h). In the IO sector, the area average wind stress during the LGM was 5 to 6% stronger than during the H1, YD, and O_H ( Figure 4). In the PO sector, the area average wind stress in the LGM was 4%, 2%, and 6% stronger than during the H1, YD, and O_H, respectively (Figure 4b-d). Table 1 shows an overview of the simulated wind stresses. The simulated upwelling core overlapped spatially with the wind stress maxima. The SO upwelling in the IO sector had a strong correlation with wind stress maxima (46°S) (r ≈ 0.8; p < 0.05) (Figure 1c,g), whereas, in the PO sector, it had a Pearson correlation coefficient of about 0.2 (p < 0.05) (Figure 1d,h). In the IO sector, the area average wind stress during the LGM was 5 to 6% stronger than during the H1, YD, and O_H ( Figure 4). In the PO sector, the area average wind stress in the LGM was 4%, 2%, and 6% stronger than during the H1, YD, and O_H, respectively (Figure 4b-d). Table 1 shows an overview of the simulated wind stresses.

Antarctic Sea Ice
The TraCE-21ka-modeled Antarctic sea ice coverage in the LGM was at a maximum in the IO sector (greater than a 70% sea ice fraction) in comparison to the PO sector (≈30% sea ice fraction) (Figure 4a). Figure 4 illustrates the deglacial decrease in the Antarctic sea

Antarctic Sea Ice
The TraCE-21ka-modeled Antarctic sea ice coverage in the LGM was at a maximum in the IO sector (greater than a 70% sea ice fraction) in comparison to the PO sector (≈30% sea ice fraction) (Figure 4a). Figure 4 illustrates the deglacial decrease in the Antarctic sea ice coverage, where its magnitude and range shrunk from the LGM to the O_H. In particular, during the O_H, the PO sector was entirely devoid of sea ice. The model simulation suggests that the SO Antarctic sea ice coverage followed a similar trend to the deglacial decrease of the Southern Hemisphere surface westerlies and wind stress. The sea ice boundary (defined as the ocean surface area's margin that is covered with more than 5% sea ice fraction) during the LGM roughly overlapped with surface density ranging from 27.5 to 27.7 kg m −3 . Figure 6 highlights the close overlap between the transition from a negative to a positive buoyancy flux (solid blue contour line represents the boundary) and the edge of the ocean surface area covered by more than an 80% sea ice fraction. The simulation showed that in both of the ocean sectors, the latitudinal position of the transition zone of the buoyancy flux shifted poleward by ≈11 • from the LGM to the O_H (Table 2). ice coverage, where its magnitude and range shrunk from the LGM to the O_H. In particular, during the O_H, the PO sector was entirely devoid of sea ice. The model simulation suggests that the SO Antarctic sea ice coverage followed a similar trend to the deglacial decrease of the Southern Hemisphere surface westerlies and wind stress. Figure 5 illustrates the equatorward expansion of the denser surface water in the LGM. The sea ice boundary (defined as the ocean surface area's margin that is covered with more than 5% sea ice fraction) during the LGM roughly overlapped with surface density ranging from 27.5 to 27.7 kg m -3 . Figure 6 highlights the close overlap between the transition from a negative to a positive buoyancy flux (solid blue contour line represents the boundary) and the edge of the ocean surface area covered by more than an 80% sea ice fraction. The simulation showed that in both of the ocean sectors, the latitudinal position of the transition zone of the buoyancy flux shifted poleward by ≈11° from the LGM to the O_H (Table 2).  ice coverage, where its magnitude and range shrunk from the LGM to the O_H. In particular, during the O_H, the PO sector was entirely devoid of sea ice. The model simulation suggests that the SO Antarctic sea ice coverage followed a similar trend to the deglacial decrease of the Southern Hemisphere surface westerlies and wind stress. Figure 5 illustrates the equatorward expansion of the denser surface water in the LGM. The sea ice boundary (defined as the ocean surface area's margin that is covered with more than 5% sea ice fraction) during the LGM roughly overlapped with surface density ranging from 27.5 to 27.7 kg m -3 . Figure 6 highlights the close overlap between the transition from a negative to a positive buoyancy flux (solid blue contour line represents the boundary) and the edge of the ocean surface area covered by more than an 80% sea ice fraction. The simulation showed that in both of the ocean sectors, the latitudinal position of the transition zone of the buoyancy flux shifted poleward by ≈11° from the LGM to the O_H ( Table 2).

Discussion
The scarcity of paleo data and differing modeling results for the SO lead to disagreements regarding the mean condition of the LGM westerly flows. A recent study used TraCE-21ka simulations to show an equatorward migration of the mean position of Southern Hemisphere westerly winds in the LGM [45], in agreement with previous studies [3]. They also showed that TraCE-21ka simulates an abrupt poleward shift during the H1 and the YD, which is consistent with the modeling and paleo proxy evidence [9,[58][59][60][61][62][63]. Moreover, in Section 2.2, we further found that the TraCE-21ka model performed well when simulating the changes in the near-surface winds and sea surface temperatures over the SO with present-day reanalysis. Therefore, we can agree that TraCE-21ka can successfully reproduce the Southern Hemisphere's natural variations during the last deglaciation.

Role of Southern Hemisphere Westerlies
TraCE-21ka simulated a pronounced zonal asymmetry during the LGM with an austral winter double jet at the SO's latitudes, which agrees with earlier studies [26,64]. This asymmetry in the Southern Hemisphere's westerly winds weakened during the H1 and the YD and was replaced by a single mid-latitude jet during the onset of the Holocene. Previous studies have mostly focused on the deglacial migration of the Southern Hemisphere's subtropical westerly jet and proxy moisture data as an analog to show the deglacial change of the zonal average Southern Hemisphere westerlies' circulation. However, it is essential to study the change in the strength of the surface Southern Hemisphere westerlies instead of the general Southern Hemisphere westerlies' circulation. Moreover, it is crucial to move beyond the average zonal perspective of the SO circulation.
Former modeling studies have shown that surface zonal wind stress changes represent the simulated change in the Southern Hemisphere westerlies' circulation [65]. Studies have revealed an increase or poleward shift of wind stress forcing results in enhanced winddriven SO upwelling [3][4][5]7,9,61,62,[66][67][68][69][70]. The simulated deglacial Southern Hemisphere surface westerly winds varied in strength and spatial disposition, imposing varying wind stress over the SO and altering the wind-driven SO upwelling. Therefore, it is necessary to study the deglacial Southern Hemisphere surface westerly winds' change over the SO.
Our result shows a strong correlation between the SO upwelling and wind stress maxima changing through the last deglaciation. In the IO sector, it decreased from very strong (r ≈ 0.8; p < 0.05) (Figure 1g) during the LGM and early deglaciation (20 to 14.9 kyr BP) to strong (r ≈ 0.7; p < 0.05) (Figure 1c) during the late deglaciation (12.8 to 9 kyr BP). However, in the PO sector, it decreased from very strong (r ≈ 0.9; p < 0.05) (Figure 1h) during the LGM and early deglaciation to strongly negative (r ≈ −0.7; p < 0.05) (Figure 1d) during the late deglaciation. Thus, our analysis supports the idea that the Southern Hemisphere surface westerlies primarily drove the SO upwelling. The simulated Southern Hemisphere surface westerlies were the weakest during the onset of the Holocene, yet the SO upwelling increased in both sectors. Therefore, the wind-driven upwelling failed to explain the anomalous upwelling during the onset of Holocene, suggesting the role of other factors in determining the deglacial change in upwelling.

Role of Sea Ice
Production of sea ice affects the seawater's density and buoyancy, and accordingly, the SO stratification and circulation [19,39,81]. Sea ice formation results in salt flux injection (wintertime brine rejection) southward and freshwater flux (summertime melting) northward of the sea ice boundary. Less buoyant, denser ocean surface water forms via brine rejection and sinks around Antarctica under the sea ice. In this region, the heat flux is weak, and instead of temperature, salt controls the ocean water density [82]. The brine rejections and sea ice melting are the biggest contributors to surface freshwater forcing. It controls the change in buoyancy (b = −g(ρ − ρ o )), which is a deviation of ocean density (ρ) from its reference value (ρ o = 1027 kg m −3 ) times gravitational acceleration (g).
Buoyancy flux (B o ) is the density flux across the air-sea interface and can be expressed as [83]: where g is the gravitational acceleration, ρ is ocean reference density, c p is the specific heat for seawater, α is the effective thermal expansion coefficient −ρ −1 ∂ρ ∂T , T is the sea surface temperature, Q o is the net surface heat flux (incoming shortwave solar radiation + outgoing sensible heat flux + latent heat flux + longwave radiation), β is the effective haline contraction coefficient ρ −1 ∂ρ ∂S , S o is the ocean surface salinity, E is the rate of evaporation, and P is the rate of precipitation. A positive buoyancy flux obtained from the definition, as shown in Equation (1) from [83] (represented by the shaded region in Figures 5 and 6) implies less buoyant surface water (buoyancy loss, i.e., densification), which is generated either by surface cooling, evaporation, or brine rejection that is associated with sea ice formation.
The 27.6 kg m −3 ocean density isopycnals coarsely mark the average transition between the lower and upper meridional overturning cells within the SO circulation [20]. Our analysis showed that the simulated deglacial sea ice boundary coarsely coincided with the surface ocean density ranging between 27.5 and 27.7 kg m −3 , as highlighted in Figure 5. Therefore, we reasoned that the simulated deglacial sea ice boundary (margin of the ocean surface area covered with more than a 5% sea ice fraction) approximately marked the boundary between the SO upper and lower meridional overturning cells'. Our investigation indicated that the edge of the ocean surface area covered with more than an 80% sea ice fraction thoroughly overlapped with the transition between the negative and positive buoyancy flux zones, as represented by a solid blue contour line in Figure 6. This implies that denser surface ocean water lay poleward to the ocean surface area covered with more than an 80% sea ice fraction, as illustrated in Figure 5.
However, the transition between the negative and positive buoyancy flux zones in the present-day climate coincided with the ocean surface area covered with a 70% sea ice fraction, as shown by [38]. This discrepancy in the simulated sea ice fraction could be due to a coarser resolution of the CSIM or the annual average data availability. Regardless, it is not significant whether the transition between the negative and positive buoyancy flux zones overlapped with 70% or 80% sea ice coverages. Our study highlighted that the transition between the negative and positive buoyancy flux zones coincided with the ocean surface covered with sea ice most of the year.
Past studies have shown that an increased salt flux during glacial times increased the Antarctic Bottom Water production and strengthened the lower overturning cell [2,38,84]. Idealized circulation model experiments [85] and dynamical scaling arguments [38] during the LGM have shown that a change in buoyancy forcing via an extended ocean sea ice coverage could expand the range of the Antarctic Bottom Water. We argue that the proportion of SO upwelling that upwells poleward to the ocean surface area covered with more than an 80% sea ice fraction would lose buoyancy by gaining surface water density, move poleward, and sink to form deep waters. Therefore, in an equatorward shift of annual ocean sea ice coverage, a larger fraction would upwell into denser surface waters, lose buoyancy, and strengthen the formation of deep waters.
In addition to an increase in deep water and Antarctic Bottom Water production, salt flux transport would also increase into the deep ocean. The passage of excess salt would make deep ocean water saltier. Moreover, during the glacial periods, extensive sea ice production in the SO would steepen the surface meridional (north-south) and the vertical (depth) density gradient [39]. The deep ocean is cold and the increase in salt flux during glacial times would stratify the ocean density structure and increase the vertical salinity gradient. The increase in a vertical salinity gradient would inhibit the vertical mixing of deep ocean waters. Furthermore, a freshwater flux from the summer sea ice melt would stratify the surface ocean [12,78,86]. Freshwater flux north of the permanent sea ice margin would likely move the upwelled water northward by increasing its buoyancy. The increase in the northward transport of more buoyant lighter surface waters would increase the meridional salinity gradient, thereby strengthening the upper meridional overturning cell [19].
Deglacial sea ice and surface density water variability would likely impact the energy exchange and the SO water mass production and circulation [42,43]. During an extended glacial sea ice coverage, the southward-flowing denser ocean water would circulate under sea ice. Therefore, sea ice would insulate the upper ocean from atmospheric forcing and act as a barrier to air-sea transfer by reducing the atmospheric exposure time of freshly upwelled waters [39][40][41].

Last Glacial Maximum
The TraCE-21ka model simulation showed that the deglacial Southern Hemisphere Pacific basin underwent significant atmospheric and oceanographic changes. During the LGM, the model simulated the strongest deglacial near-surface zonal winds at the SO's latitudes. Recent studies support the TraCE-21ka simulation, where they agree with the stronger LGM Southern Hemisphere surface westerlies [39,80]. The TraCE-21ka model simulated the most vigorous deglacial SO upwelling during the LGM, and a study by [87] supports the simulation. Our analysis indicated that surface westerlies primarily drove the SO upwelling.
During the LGM, global cooling due to lower atmospheric greenhouse gas concentrations (Figure 1e) resulted in the simulated equatorward expansion of the glacial sea ice coverage area to ≈45 • S in the IO and PO sectors (Figure 4a). The model-simulated equatorward extension of the LGM sea ice agrees with previous studies [38,79,80,86,[88][89][90]. Our results highlighted that during the LGM, the ocean surface area covered with more than an 80% sea ice fraction overlapped with a transition between the negative and positive buoyancy fluxes ( Figure 6). It was positioned ≈11 • equatorward compared to during the O_H ( Table 2). The simulated equatorward extension of denser ocean surface water during the LGM agrees with past studies [12,38]. The simulation showed that the LGM deep ocean was extensively stratified (Figure 7a) and colder and saltier (Figure 7b) than during the O_H and agreed with previous studies [12,28,[91][92][93].
The simulation showed that the extensive coverage of sea ice during the LGM led to a massive influx of salt into the ocean compared to the O_H, as shown in Figure 8a. Our results showed that the vertical salinity (illustrated in Figure 7a) and density gradient of the LGM deep ocean water was seven and two times stronger than during the O_H, respectively, as shown in Table 3. This extensive stratification of the LGM deep ocean waters would subsequently lead to low vertical diffusion and less mixing of the LGM waters [92,94]. Regarding the LGM, Figure 8c highlights the intensive freshwater flux into lighter surface waters at the sea ice boundary. Our analysis showed that the LGM meridional salinity gradient was stronger by a factor of nine relative to during the O_H (Table 3). This increase in simulated freshwater flux and meridional salinity gradient would strengthen the glacial upper meridional overturning circulation, which is consistent with [19]. The simulation of the Global Meridional Overturning Circulation (0.5 to 4.5 km depth) during the LGM (Figure 9a) showed an extended lower meridional overturning cell (blue), an upper meridional overturning cell (red), and a shallower Deacon cell (≈45 • S to 60 • S; red). The LGM lower meridional overturning cell comprising the Antarctic Bottom Water and its derivatives expanded above the average topographic features in the SO to depths of about 1600 m (Figure 9b). Previous studies have also reported expanded lower meridional overturning circulation during the LGM [38,90,95,96]. Our analysis showed that the strength of the maximum lower meridional overturning transport during the LGM was about 1.3 times stronger than during the O_H (Table 3).  The simulation showed that the extensive coverage of sea ice during the LGM led to a massive influx of salt into the ocean compared to the O_H, as shown in Figure 8a. Our results showed that the vertical salinity (illustrated in Figure 7a) and density gradient of the LGM deep ocean water was seven and two times stronger than during the O_H, respectively, as shown in Table 3. This extensive stratification of the LGM deep ocean waters would subsequently lead to low vertical diffusion and less mixing of the LGM waters [92,94]. Regarding the LGM, Figure 8c highlights the intensive freshwater flux into lighter surface waters at the sea ice boundary. Our analysis showed that the LGM meridional salinity gradient was stronger by a factor of nine relative to during the O_H (Table 3). This increase in simulated freshwater flux and meridional salinity gradient would strengthen the glacial upper meridional overturning circulation, which is consistent with [19]. The simulation of the Global Meridional Overturning Circulation (0.5 to 4.5 km depth) during the LGM (Figure 9a) showed an extended lower meridional overturning cell (blue), an upper meridional overturning cell (red), and a shallower Deacon cell (≈45°S to 60°S; red). The LGM lower meridional overturning cell comprising the Antarctic Bottom Water and its derivatives expanded above the average topographic features in the SO to depths of about 1600 m (Figure 9b). Previous studies have also reported expanded lower meridional overturning circulation during the LGM [38,90,95,96]. Our analysis showed that the strength of the maximum lower meridional overturning transport during the LGM was about 1.3 times stronger than during the O_H (Table 3).  Note: MSG-meridional salinity gradient, VSG-vertical salinity gradient, VDG-vertical density gradient, IO-Indian Ocean, PO-Pacific Ocean, LGM-Last Glacial Maximum, O_H-onset of the Holocene, MOC-Meridional Overturning Circulation. The meridional salinity gradient is the difference between the salinity within 65 • S to 54 • S and 54 • S to 43 • S above a 70-meter ocean depth. The vertical gradients are the differences between 1.8 to 2.7 km and 2.7 to 4.7 km ocean depths within 50 • S to 30 • S. The salinity gradient unit is in grams per kilogram (g kg −1 ), and the density gradient is in kilograms per cubic meter (kg m −3 ). The maximum lower meridional overturning transport is the model maximum meridional transport below 2 km of ocean depth; its unit is in Sverdrups (Sv).  Note: MSG-meridional salinity gradient, VSG-vertical salinity gradient, VDG-vertical density gradient, IO-Indian Ocean, PO-Pacific Ocean, LGM-Last Glacial Maximum, O_H-onset of the Holocene, MOC-Meridional Overturning Circulation. The meridional salinity gradient is the difference between the salinity within 65°S to 54°S and 54°S to 43°S above a 70-meter ocean depth. The vertical gradients are the differences between 1.8 to 2.7 km and 2.7 to 4.7 km ocean depths within 50°S to 30°S. The salinity gradient unit is in grams per kilogram (g kg -1 ), and the density gradient is in kilograms per cubic meter (kg m -3 ). The maximum lower meridional overturning transport is the model maximum meridional transport below 2 km of ocean depth; its unit is in Sverdrups (Sv). The expansion of the LGM sea ice resulted in the extension and shoaling of the lower meridional overturning cell, as shown in Figure 9b. The deep ocean waters were saltier and extensively stratified. This led to weaker vertical diffusion and a lack of turbulent mixing of carbon-rich waters from the lower meridional overturning cell with the overlying carbondeficient waters [38]. The poorly mixed water would reach the ocean surface under the extended coverage of sea ice, which would likely restrict the outgassing of CO 2 [42,43,97]. Furthermore, the TraCE-21ka-simulated colder and expanded LGM deep ocean water would increase the deep ocean carbon storage [96].

Deglaciation: H1 and YD
The wind stress-upwelling analogy was retained well in the TraCE-21ka-simulated H1 and YD climatologies. During the H1 and the YD, the model simulated a southward shift in the Intertropical Convergence Zone position ( Figure S3) and more vigorous wind stress forcing in the latitudes of the SO in comparison to the O_H. The simulated increase in the SO upwelling in response to the strengthened surface Southern Hemisphere westerlies during the H1 and the YD was consistent with the upwelling records [9,98].
The deglacial climate experienced Northern Hemisphere meltwater discharge during the H1 and YD (Figure 1b,f) events. The simulation showed that the northern Atlantic meltwater discharge weakened the Atlantic Meridional Ocean Circulation and reduced the supply of North Atlantic origin waters into the SO's latitudes. The meltwater forcings led to the associated Southern Hemisphere Pacific bipolar seesaw [99], and deep ocean warming resulted in the Antarctic sea ice's deglacial retreat. During the H1 and the YD, the sea ice retreated poleward from its equatorward position during the LGM. The SO deep waters were fresher and less stratified than the LGM waters. The simulation showed that the lower overturning cell during the H1 and the YD was less expanded than the LGM. Consequently, this may have resulted in more turbulent mixing of carbon-rich water from the lower meridional overturning cell with the overlying waters [12,38]. The expansion of the LGM sea ice resulted in the extension and shoaling of the lower meridional overturning cell, as shown in Figure 9b. The deep ocean waters were saltier and extensively stratified. This led to weaker vertical diffusion and a lack of turbulent mixing of carbon-rich waters from the lower meridional overturning cell with the overlying carbon-deficient waters [38]. The poorly mixed water would reach the ocean surface under the extended coverage of sea ice, which would likely restrict the outgassing of CO2 [42,43,97]. Furthermore, the TraCE-21ka-simulated colder and expanded LGM deep ocean water would increase the deep ocean carbon storage [96] The wind stress-upwelling analogy was retained well in the TraCE-21ka-simulated H1 and YD climatologies. During the H1 and the YD, the model simulated a southward shift in the Intertropical Convergence Zone position ( Figure S3) and more vigorous wind stress forcing in the latitudes of the SO in comparison to the O_H. The simulated increase Enhanced upwelling during the H1 and the YD would likely bring more carbonrich water to the ocean's surface. The reduction in the sea ice coverage would enhance the outgassing of CO 2 [42,43,97]. Furthermore, the ocean water during the H1 and the YD was warmer than the LGM and decreased the solubility of CO 2 . The TraCE-21kamodeled upwelling records during the H1 and the YD were consistent with the increase in atmospheric CO 2 concentrations from ice core records [11] (Figure 1a,e).

Deglaciation: Onset of the Holocene
The wind stress-upwelling analogy did not fit with the inconsistent increase in upwelling during the O_H, especially in the PO sector (Figure 1d). Meltwater fluxes from changing continental ice sheets and increases in insolation and greenhouse gasses resulted in Southern Hemisphere deglacial deep ocean warming (Figure 7b) and Antarctic sea ice's deglacial retreat (Figure 4). The simulated IO sector was partially devoid of sea ice and the PO sector was entirely free of sea ice (Figure 4d) during the O_H.
Our results highlighted that the absence of sea ice made the SO deep waters fresher ( Figure 7b) and less stratified (Figure 7a) compared to the LGM. The lower meridional overturning cell's upper margin during the O_H moved deeper to about 1900 m than the shallower LGM lower meridional overturning cell (Figure 9b). The simulation showed that the strength of the upper meridional overturning cell (Atlantic Meridional Ocean Circulation) was about two times stronger during the O_H compared with during the H1 and the YD. This would likely bring low carbon-rich waters into the upper meridional overturning cell. Additionally, the magnitude of the simulated freshwater flux increased by 3 times in the IO and 52 times in the PO sectors compared to the LGM, as shown in Figure 8d. The increased freshwater flux at the O_H sea ice boundary latitudes would strengthen the upper meridional overturning cell's circulation. Therefore, our analysis highlighted that the increase in freshwater flux during the O_H explains the anomalous SO upwelling.
Modern satellite and modeling studies have emphasized the importance of the seaice-induced freshwater flux. Recent studies have indicated that the SO freshwater flux alters the salinity, vertical and meridional density gradients, and thermohaline circulation of the modern ocean [100]. Moreover, an estimated annual northward sea ice freshwater flux from Antarctica [100] shows comparable sea ice freshwater transport to the O_H. Our study suggested that it is vital to understand the implication of sea ice induced freshwater transport changes to understand the future SO climate changes.
Moreover, the absence of sea ice would increase energy transfer efficiency from the atmosphere to the ocean [42,43] compared to the LGM, suggesting an increase in SO upwelling. The dynamical changes in the carbon content of SO upwelling water are beyond the scope of this study. However, it is crucial to understand that the deglacial retreat of Antarctic sea ice changed the SO circulation and the vertical and turbulent mixing of SO deep waters, thereby changing the upwelling waters' carbon content.

Conclusions
This article investigated the deglacial Southern Ocean upwelling's evolution on a subbasin scale using a fully coupled global climate TraCE-21ka model. The main conclusions from this work are as follows: (1). The model simulated the deglacial evolution of the Southern Hemisphere surface westerlies, and the Antarctic sea ice, which agreed with proxy records. It established TraCE-21ka's consistency and applicability in simulating the physical mechanisms that determine the Southern Ocean upwelling. (2). The Southern Hemisphere surface westerly winds were strong during the Last Glacial Maximum, Heinrich Stadial 1, and Younger Dryas compared to during the onset of the Holocene, with the Last Glacial Maximum's westerly winds being the strongest. (3). Intense surface westerly winds and the subsequent wind stress forcing over the Southern Ocean resulted in more vigorous upwelling during the Last Glacial Maximum, Heinrich Stadial 1, and Younger Dryas. Our investigation indicated that the strength of the wind stresses primarily governs the variability in the Southern Ocean upwelling around Antarctica. (4). The strong simulated upwelling in both the ocean sectors during the onset of the Holocene cannot be explained by weaker wind stress forcing. This study suggested that the sea ice modulated the buoyancy forcing via an enhanced freshwater flux, and together with the zonal wind stress forcing, describes the strong upwelling during the onset of Holocene. (5). The Antarctic sea ice expansion during the Last Glacial Maximum was dynamically linked to an increase in deep ocean stratification. Moreover, it resulted in the extension, shoaling, and increase in transport of the lower meridional overturning circulation (Antarctic Bottom Water and its derivatives) during the Last Glacial Maximum. Additionally, the Antarctic sea ice affected the surface buoyancy flux via freshwater and salt fluxes. It also regulated the Southern Ocean circulation by marking the division between the upper (Atlantic Meridional Ocean Circulation) and lower meridional overturning ocean circulation. (6). This study highlighted that the evolution of the Southern Hemisphere westerly winds and sea ice history is essential in understanding the ocean-atmosphere coupling at high latitudes and its role in changing the deglacial Southern Ocean upwelling. The TraCE-21ka simulation demonstrated that the Southern Ocean atmosphere, ocean, and cryosphere were dynamically interlinked. The changes in the Southern Hemisphere westerlies, Southern Ocean upwelling, and Antarctic sea ice change must have cooccurred through the last deglaciation during the late Pleistocene. This relationship is vital to include in future climate modeling studies to understand the Southern Ocean upwelling variations. Our results may help to understand the natural global warming process and its consequences for the Southern Ocean dynamics.