The Roles of Orbital and Meltwater Climate Forcings on the Southern Ocean Dynamics during the Last Deglaciation

: The last deglacial climate evolution, from 19 to 9 thousand years before the present, represents the vital role of feedback in the Earth’s climate system. The Southern Ocean played a fundamental role by exchanging nutrients and carbon-rich deep ocean water with the surface during the last deglaciation. This study employs a fully coupled Earth system model to investigate the evolution of Southern Ocean dynamics and the roles of changes in orbital and meltwater forcings during the last deglaciation. The simulation supports that the Southern Ocean upwelling was primarily driven by windstress. The results show that the melting and formation of Antarctic sea ice feedback inﬂuenced Southern Ocean surface buoyancy ﬂux. The increase in Antarctic sea ice melt-induced freshwater ﬂux resulted in a steepened north-south surface salinity gradient in the Southern Ocean, which enhanced the upwelling. The single-forcing experiments indicate that the deglacial changes in orbital insolation inﬂuenced the Southern Ocean upwelling. The experiments also highlight the dominant role of Northern Hemisphere meltwater discharge in the upper and lower branch of the Meridional Overturning Circulation. Furthermore, orbital forcing shows lesser deglacial Antarctic sea ice retreat than the Northern Hemisphere meltwater forcing, which follows the bipolar seesaw mechanism.


Introduction
Numerous studies have shown that changes in orbital insolation played a significant role in driving the Earth's glacial-interglacial cycles [1][2][3]. For example, previous studies have suggested that an increase in mid-latitude to high-latitude Northern Hemisphere (NH) spring-summer insolation triggered the last deglaciation [3][4][5][6]. Consequently, the meltwater discharge from continental ice sheets retreat weakened the Atlantic Meridional Overturning Circulation (AMOC) [6] and caused the Southern Hemisphere (SH) deglacial warming via the bipolar seesaw mechanism [7][8][9][10]. The orbitally driven changes in NH insolation are also associated with increased Antarctic temperatures and atmospheric carbon dioxide concentrations during the last deglaciation [11]. Therefore, it is vital to understand the role of orbital forcing and its feedback in order to understand the Earth's climate.
The ocean has played a fundamental role in the Earth's climate system by absorbing, transporting, and releasing heat and carbon. The Southern Ocean (SO) (south of 30 • S) ventilates the deep ocean by exchanging heat and gases between the deep ocean and the atmosphere [12]. The deglacial period from 19 to 9 kyr (thousand years) BP (before present) is punctuated by abrupt changes in SO ventilation [13][14][15][16], with about 75 parts per million by volume increase in atmospheric carbon dioxide concentrations [17,18]. Therefore, the SO has played a crucial role in controlling past atmospheric carbon dioxide concentrations [19]. SO climate during the last deglaciation. Recent studies have also highlighted the role of orbital forcing in wind-driven SO upwelling [76,77]. This modeling study investigates the response of orbital and meltwater climate forcings to the SO dynamics during the last deglacial period from 19 to 9 kyr BP. Previous studies have indicated that local SO upwelling exhibits distinctive evolution patterns than the zonal mean SO upwelling during the last deglaciation [58]. This study examines the evolution of SO upwelling in the Atlantic Ocean basin to understand the MOC in the SO and North Atlantic. This contribution also focuses on the evolution of SO upwelling in response to Antarctic sea ice melt-induced freshwater fluxes using a fully coupled Earth system model during the last deglaciation. The conclusions support that the SO upwelling is influenced by orbital forcing and is primarily driven by SH windstress. The increase in Antarctic sea ice melt-induced freshwater export on the surface of the SO enhanced the SO upwelling during the last deglacial period.

Data
We analyzed annual data from the TraCE-21ka climate model for the last deglacial period. The TraCE-21ka is a fully coupled land surface-atmosphere (Community Atmospheric Model 3)-ocean (Parallel Ocean Program)-sea ice (Community Sea Ice Model) climate model experiment using the Community Climate System Model version 3 [78]. The TraCE-21ka ice and ocean data have an identical resolution (gx3v5; about 3.6 degrees longitudinal and a variable latitudinal resolution to about 0.9 degrees near the equator). The ice and ocean data were interpolated to ≈3.75 degrees resolution identical with the atmospheric data (T31). The ocean model has a vertical z coordinate with 25 depth levels and uses the ocean eddies parameterization [79]. In addition, the dynamic thermodynamic sea ice model incorporates a subgrid-scale ice thickness distribution.
We analyzed complete (FULL), orbital (ORB), and meltwater (MWF) single-forcing transient simulations for the last deglacial period. The focus was to investigate FULL, ORB, and MWF climate forcings' contribution to the upwelling and MOC in the SO. The experimental data analyzed are: (1) FULL: a complete model simulation of the past 22 kyr BP with transient alteration in orbital insolation, meltwater fluxes, continental ice sheets, and greenhouse gases. The continental ice sheet extents and orographies were adopted from the ICE-5G (VM2) [80].  Figure 1b) [81]. The model also simulated the meltwater pulse 1A (mwp-1A) at about 14.1 kyr BP, which caused the Allerød warm period (≈14 to 12.9 kyr BP). The mwp-1A is described by the sea level rise of about 20 m in about 300 to 500 years ( Figure 1b) [82]. The author of [80] assigned the total 20 m sea level rise to NH origin, which set off an AMOC collapse in climate models [83]. However, several studies have included meltwater forcing contributions from the Antarctic (SH) (15 m of equivalent sea level volume) and Laurentide Ice Sheets (NH) (5 m of equivalent sea level volume) [84,85]. Studies have also shown that a NH meltwater forcing of more than 5 m resulted in a complete shutdown of AMOC, which is inconsistent with proxy records [86]. Therefore, the TraCE-21ka adopted the SH (Antarctic) meltwater forcing (60 Sverdrup (Sv); 1 Sv = 10 6 m 3 s −1 ) and the NH meltwater forcing (20 Sv) during the mwp-1A. This meltwater discharge was added as a freshwater flux to the ocean model onto the ocean surface.

Model Performance
Numerous studies of the last deglacial period have shown that the TraCE-21ka model reproduced several critical features of the global climate evolution [78,87]. For example, the TraCE-21ka model replicated the Asian-African monsoon [88], El Niño [89], the highlatitude seasonal temperature [90], the AMOC transport, and the SH regional sea surface temperature [5]. In addition, the model reasonably simulated the present-day atmospheric circulation and sea surface temperature in the SO [49,91].
We performed a comprehensive data-model comparison between proxy records (Figure 1d,e) and the TraCE-21ka FULL simulation upwelling in the Atlantic Ocean sector of the SO during the last deglacial period (Figure 1e). The TraCE-21ka model simulated an increase in SO upwelling during the Heinrich 1 (H1) and the latter part of the Younger Dryas (YD) events in agreement with numerical simulations [15] and proxy records (Figure 1d,e) [13,14,18,92]. The simulation demonstrated a gradual increase in SO upwelling during the early part of deglaciation, agreeing with the proxy record [13]. As such, it did not replicate the abrupt increase around 16.3 kyr BP, as shown by proxy record [14] and numerical simulations [15]. We found that the simulated SO upwelling followed the high precision, century-scale proxy record from [14] reasonably well ( Figure 1e). The model simulation considerably overlapped with the century-scale changes in δ 13 C − CO 2 proxy marine source records at 12.9 kyr BP, indicating strengthened air-sea gas exchange in the SO [14]. The simulation also showed an increase in the SO upwelling around 17 to 16 kyr BP associated with increased atmospheric carbon dioxide concentrations, evident from the proxy record [18] and numerical simulations [15]. Additionally, the authors of [58] employed the TraCE-21ka model and simulated SO upwelling at corresponding grids of proxy cores in the SO [13,93]. Their simulation also followed the proxy records, with some alterations in the magnitude and timing of the deglacial upwelling changes. In conclusion, the TraCE-21ka model performed reasonably well in simulating the changes in the SO upwelling during the last deglacial period.
Like most global ocean models, the TraCE-21ka model did not simulate the eddy effect. Instead, the ocean model used the Gent and McWilliams eddy parameterization [79] to approximate the eddy-induced vertical velocity response to surface forcings. Thus, the TraCE-21ka model produced an Eulerian-mean meridional overturning circulation. Nevertheless, it performed reasonably well in representing the deglacial SO upwelling due to a significant contribution from Eulerian-mean upwelling to tracer movements.

Methodology
The authors of [58] employed the TraCE-21ka model and found that the zonal mean upwelling from vertical ocean velocity equals the zonal mean Ekman pumping in the upper 2000 m in the SO. This study diagnosed the TraCE-21ka data and integrated the upper 2000 m of vertical ocean velocity to calculate the SO upwelling volume in the Atlantic Ocean basin during the deglacial period, as shown in Equation (1).
whereŴ is the integrated upwelling along with the depth, w(k) is the vertical velocity at each ocean grid, dz(k) is the thickness of the kth level, and N is the number of vertical ocean grids. We also analyzed the integrated upwelling for the upper 500 m of the SO. Similar conclusions were drawn from the study despite a change in depth of the representative sector. This study examined the time evolution of the SO upwelling through the last deglacial period from 19 to 9 kyr BP. It included the 19 kyr BP, the H1, and the YD events. We selected two millennial timescales: (1) the Heinrich 1 (H1; from 17.3 to 16.3 kyr BP), and (2) the onset of the Holocene (O_H; from 10 to 9 kyr BP), to highlight early and later parts of the deglaciation.

Southern Ocean Dynamics
We analyzed the TraCE-21ka complete (FULL) and single forcing orbital (ORB) and meltwater (MWF) simulations during the last deglaciation to understand the evolution of the MOC. Figure 2a illustrates that the TraCE-21ka FULL experiment indicated a weakened upper overturning cell during the early from the later part of deglaciation. It also shows a more robust and extended lower overturning cell during the early part compared to the later part of the deglaciation (see Supplementary Figure S1). We found that only the MWF experiment shows weakened AMOC similar to the FULL experiment during the early part of the deglaciation (Figure 2b,c).

Methodology
The authors of [58] employed the TraCE-21ka model and found that the zonal mean upwelling from vertical ocean velocity equals the zonal mean Ekman pumping in the upper 2000 m in the SO. This study diagnosed the TraCE-21ka data and integrated the upper 2000 m of vertical ocean velocity to calculate the SO upwelling volume in the Atlantic Ocean basin during the deglacial period, as shown in Equation (1).
where Ŵ is the integrated upwelling along with the depth, w(k) is the vertical velocity at each ocean grid, dz(k) is the thickness of the kth level, and N is the number of vertical ocean grids. We also analyzed the integrated upwelling for the upper 500 m of the SO. Similar conclusions were drawn from the study despite a change in depth of the representative sector. This study examined the time evolution of the SO upwelling through the last deglacial period from 19 to 9 kyr BP. It included the 19 kyr BP, the H1, and the YD events. We selected two millennial timescales: (1) the Heinrich 1 (H1; from 17.3 to 16.3 kyr BP), and (2) the onset of the Holocene (O_H; from 10 to 9 kyr BP), to highlight early and later parts of the deglaciation.

Southern Ocean Dynamics
We analyzed the TraCE-21ka complete (FULL) and single forcing orbital (ORB) and meltwater (MWF) simulations during the last deglaciation to understand the evolution of the MOC. Figure 2a illustrates that the TraCE-21ka FULL experiment indicated a weakened upper overturning cell during the early from the later part of deglaciation. It also shows a more robust and extended lower overturning cell during the early part compared to the later part of the deglaciation (see Supplementary Figure S1). We found that only the MWF experiment shows weakened AMOC similar to the FULL experiment during the early part of the deglaciation (Figure 2b,c).   The FULL experiment reasonably simulated the decrease in AMOC transport during the last deglaciation, which is consistent with the AMOC export proxy record (Figure 3a). We found that the simulated strength of upper overturning circulation during the H1 was about 57% weaker than during the O_H (Table 1). Conversely, the simulated maximum strength of the lower overturning circulation during the H1 was about 12% stronger than it was during the O_H (Figure 3b) ( Table 1).
The FULL experiment reasonably simulated the decrease in AMOC transport during the last deglaciation, which is consistent with the AMOC export proxy record (Figure 3a). We found that the simulated strength of upper overturning circulation during the H1 was about 57% weaker than during the O_H (Table 1). Conversely, the simulated maximum strength of the lower overturning circulation during the H1 was about 12% stronger than it was during the O_H (Figure 3b) ( Table 1).  We found that the MWF experiment reasonably reproduced the strength of upper and lower overturning circulation compared to the ORB experiment during the last deglaciation (Figure 3a,b). The MWF experiment simulated strength of upper ( Figure 3a) and lower (Figure 3b) overturning circulation has an excellent agreement with the FULL experiment during the early ((r = 1, p < 0.001) and (r = 0.7, p < 0.001)) and late ((r = 0.9, p < 0.001) and (r = 0.8, p < 0.001)) parts of the deglaciation, respectively. Figure 4a indicates that the ORB experiment simulated maximum SO upwelling (the scale is on the right y-axis) agrees with the FULL experiment (the scale is on the left y-axis)  We found that the MWF experiment reasonably reproduced the strength of upper and lower overturning circulation compared to the ORB experiment during the last deglaciation (Figure 3a,b). The MWF experiment simulated strength of upper ( Figure 3a) and lower (Figure 3b) overturning circulation has an excellent agreement with the FULL experiment during the early ((r = 1, p < 0.001) and (r = 0.7, p < 0.001)) and late ((r = 0.9, p < 0.001) and (r = 0.8, p < 0.001)) parts of the deglaciation, respectively. Figure 4a indicates that the ORB experiment simulated maximum SO upwelling (the scale is on the right y-axis) agrees with the FULL experiment (the scale is on the left y-axis) reasonably well during the early (decrease) and late (increase) parts of the deglaciation. Specifically, the ORB experiment has a Pearson correlation coefficient of 0.4 (p < 0.01) during the early (19 to 14.7 kyr BP) and late (14.7 to 12.4 kyr BP) parts of deglaciation. On the contrary, the MWF experiment (the scale is on the left y-axis) indicated that the SO upwelling followed the FULL experiment only during the middle (14.7 to 12.4 kyr BP) part of the deglaciation (r = 0.5; p < 0.001), as shown in Figure 4a.
reasonably well during the early (decrease) and late (increase) parts of the deglaciation. Specifically, the ORB experiment has a Pearson correlation coefficient of 0.4 (p < 0.01) during the early (19 to 14.7 kyr BP) and late (14.7 to 12.4 kyr BP) parts of deglaciation. On the contrary, the MWF experiment (the scale is on the left y-axis) indicated that the SO upwelling followed the FULL experiment only during the middle (14.7 to 12.4 kyr BP) part of the deglaciation (r = 0.5; p < 0.001), as shown in Figure 4a.
The FULL experiment demonstrated that the evolution of the SO upwelling followed the changes in SH windstress reasonably well throughout the last deglaciation (Figure 4b).
We found that the simulated SO upwelling and zonal windstress has a Pearson correlation coefficient of 0.6 (p < 0.001) during the last deglacial period (19 to 9 kyr BP). Specifically, it has a Pearson correlation coefficient of 0.9 (p = 0) during the early, 0.7 (p < 0.001) during the middle, and 0.4 (p < 0.001) during the late parts of the deglaciation. Additionally, the ORB experiment demonstrated a Pearson correlation coefficient of 0.8 (p < 0.001) between SO upwelling and zonal windstress during the last deglacial period.

Antarctic Sea Ice Dynamics
The FULL experiment simulated a decrease in Antarctic sea ice coverage during the last deglacial period (Figures 5 and S2a). The FULL and the MWF experiments demonstrated a steep decrease in Antarctic sea ice coverage during the H1, the YD events, and an abrupt increase during the Antarctic Cold Reversal (ACR; ≈14.5 to ≈12.5 kyr BP) period. On the other hand, the ORB experiment simulated an extended deglacial Antarctic sea ice The FULL experiment demonstrated that the evolution of the SO upwelling followed the changes in SH windstress reasonably well throughout the last deglaciation (Figure 4b).
We found that the simulated SO upwelling and zonal windstress has a Pearson correlation coefficient of 0.6 (p < 0.001) during the last deglacial period (19 to 9 kyr BP). Specifically, it has a Pearson correlation coefficient of 0.9 (p = 0) during the early, 0.7 (p < 0.001) during the middle, and 0.4 (p < 0.001) during the late parts of the deglaciation. Additionally, the ORB experiment demonstrated a Pearson correlation coefficient of 0.8 (p < 0.001) between SO upwelling and zonal windstress during the last deglacial period.

Antarctic Sea Ice Dynamics
The FULL experiment simulated a decrease in Antarctic sea ice coverage during the last deglacial period ( Figure 5 and Figure S2a). The FULL and the MWF experiments demonstrated a steep decrease in Antarctic sea ice coverage during the H1, the YD events, and an abrupt increase during the Antarctic Cold Reversal (ACR; ≈14.5 to ≈12.5 kyr BP) period. On the other hand, the ORB experiment simulated an extended deglacial Antarctic sea ice coverage which gradually decreased during the last deglacial period. The ORB experiment showed reduced coverage of the Antarctic sea ice in the Indian and the Pacific sectors of the SO (Supplementary Figure S2b). However, the MWF experiment showed similar Antarctic sea ice coverage between the early and late parts of the deglaciation (Supplementary Figure S2c).  Figure S2b). However, the MWF experiment showed similar Antarctic sea ice coverage between the early and late parts of the deglaciation (Supplementary Figure S2c).  Figure 6 shows the decrease in the Antarctic quasi-permanent sea ice coverage (SO surface area covered with more than eighty percent sea ice fraction) in the FULL experiment during the last deglacial period. The SO region with negative surface buoyancy flux (denser surface ocean density) was also displaced towards Antarctica during the last deglacial period. Additionally, Figure 6 highlights that the Antarctic quasi-permanent sea ice boundary overlaps with the transition of surface buoyancy flux from positive (surface buoyancy gain) to negative (surface buoyancy loss) during the last deglacial period.  Figure 6 shows the decrease in the Antarctic quasi-permanent sea ice coverage (SO surface area covered with more than eighty percent sea ice fraction) in the FULL experiment during the last deglacial period. The SO region with negative surface buoyancy flux (denser surface ocean density) was also displaced towards Antarctica during the last deglacial period. Additionally, Figure 6 highlights that the Antarctic quasi-permanent sea ice boundary overlaps with the transition of surface buoyancy flux from positive (surface buoyancy gain) to negative (surface buoyancy loss) during the last deglacial period.

Discussion
This study examined the deglacial evolution of the SO upwelling in the Atlantic Ocean basin to understand the SO dynamics employing a fully coupled climate model simulation, the TraCE-21ka. Numerous studies have suggested that the SO upwelling is primarily wind-driven [12,13,15,23,94]. The authors of [23] proposed the "westerlies shift" hypothesis, wherein a poleward (equatorward) shift or an increase (decrease) in SH westerlies increased (decreased) the wind-driven SO upwelling. As a result, it allowed more (less) respired carbon dioxide to outgas into the atmosphere. Recently, the authors of [49] showed that the TraCE-21ka model simulated the deglacial shift of the SH westerly jet in agreement with the "westerlies shift" hypothesis. Moreover, the authors of [58] also employed the TraCE-21ka model and found that the local SO upwelling is driven by surface wind stress and surface buoyancy forcing and is influenced by the local topography. This work studied the basin-scale SO upwelling in the Atlantic Ocean sector during the last deglacial period instead of local upwelling. We found that the changes in the SO upwelling reasonably agreed with the proxies and numerical simulations and are driven by surface windstress and Antarctic sea ice-influenced surface buoyancy forcing.
The Antarctic sea ice played a pivotal role in the MOC by modifying water masses globally [60]. The Antarctic sea ice-ocean feedback included wintertime sea ice formation (brine rejection), which left salt content, and included summertime melting, which returned liquid freshwater to the ocean. This feedback altered surface buoyancy fluxes and modulated SO deep ocean stratification and overturning circulation [22,24,57,58,[61][62][63][64][65][66]69,70,95]. The TraCE-21ka FULL experiment indicated a deglacial decrease in Antarctic sea ice coverage in good agreement with the proxy records [96]. This study demonstrated that, during the last deglacial period, the boundary of SO surface area covered with sea ice most of the year overlapped with the transition of surface buoyancy flux (positive (gain) to negative (loss)). The authors of [61] showed a similar relationship between the Antarctic quasi-permanent sea ice and surface buoyancy flux in the present-day and the Last Glacial Maximum climate. However, this study highlighted the association between the Antarctic sea ice coverage and the surface buoyancy forcing during the last deglacial period. Our results suggested that, during the last deglaciation, an increase (decrease) in

Discussion
This study examined the deglacial evolution of the SO upwelling in the Atlantic Ocean basin to understand the SO dynamics employing a fully coupled climate model simulation, the TraCE-21ka. Numerous studies have suggested that the SO upwelling is primarily wind-driven [12,13,15,23,94]. The authors of [23] proposed the "westerlies shift" hypothesis, wherein a poleward (equatorward) shift or an increase (decrease) in SH westerlies increased (decreased) the wind-driven SO upwelling. As a result, it allowed more (less) respired carbon dioxide to outgas into the atmosphere. Recently, the authors of [49] showed that the TraCE-21ka model simulated the deglacial shift of the SH westerly jet in agreement with the "westerlies shift" hypothesis. Moreover, the authors of [58] also employed the TraCE-21ka model and found that the local SO upwelling is driven by surface wind stress and surface buoyancy forcing and is influenced by the local topography. This work studied the basin-scale SO upwelling in the Atlantic Ocean sector during the last deglacial period instead of local upwelling. We found that the changes in the SO upwelling reasonably agreed with the proxies and numerical simulations and are driven by surface windstress and Antarctic sea ice-influenced surface buoyancy forcing.
The Antarctic sea ice played a pivotal role in the MOC by modifying water masses globally [60]. The Antarctic sea ice-ocean feedback included wintertime sea ice formation (brine rejection), which left salt content, and included summertime melting, which returned liquid freshwater to the ocean. This feedback altered surface buoyancy fluxes and modulated SO deep ocean stratification and overturning circulation [22,24,57,58,[61][62][63][64][65][66]69,70,95]. The TraCE-21ka FULL experiment indicated a deglacial decrease in Antarctic sea ice coverage in good agreement with the proxy records [96]. This study demonstrated that, during the last deglacial period, the boundary of SO surface area covered with sea ice most of the year overlapped with the transition of surface buoyancy flux (positive (gain) to negative (loss)). The authors of [61] showed a similar relationship between the Antarctic quasi-permanent sea ice and surface buoyancy flux in the present-day and the Last Glacial Maximum climate. However, this study highlighted the association between the Antarctic sea ice coverage and the surface buoyancy forcing during the last deglacial period. Our results suggested that, during the last deglaciation, an increase (decrease) in the Antarctic sea ice extent corresponds to an increase (decrease) in the coverage of the buoyancy loss zone. Recent studies have also indicated that a shift in the Antarctic sea ice extent caused a change in the extent of surface buoyancy flux [21,61,62,97].
Previous studies have indicated that the SO surface buoyancy forcing was dominated by changes in freshwater export [70] associated with summertime melting of Antarctic sea ice [22,57,69]. The authors of [57] estimated the recent (2005 to 2010) bulk freshwater volume fluxes at the SO surface south of 50 • S from ocean and ice observations. They showed that the net volume of freshwater fluxes from the Antarctic sea ice melt was almost ten times more than that of the glacial ice melt [57]. Additionally, the authors of [58] employed the TraCE-21ka model and showed that changes in surface SO freshwater fluxes were essentially contributed by Antarctic sea ice melt and brine rejection during the last deglaciation. Thus, present-day studies and last deglacial studies highlight that the freshwater discharge from Antarctic sea ice melt is vital to understanding SO dynamics.
This study analyzed the zonal freshwater flux in the SO at the boundary of the quasipermanent sea ice and buoyancy flux transition (≈50 • S). We found that the simulated freshwater discharge closely followed the Antarctic sea ice coverage in all the experiments during the last deglacial period ( Figure 5). In addition, the FULL and the MWF experiments showed maximum upwelling ( Figure 4a) and enhanced Antarctic sea ice coverage ( Figure 5) at about 14.1 kyr BP that coincided with the mwp-1A.
The FULL experiment simulated the mwp-1A, which consisted of the NH and the SH meltwater forcings (Figure 1b). The meltwater discharge was added as a freshwater flux onto the ocean surface in the ocean model. The simulated freshwater flux from the SH meltwater forcing is three times the magnitude of the NH meltwater forcing (Figure 1b). On the other hand, the MWF experiment simulated the mwp-1A with only the NH meltwater discharge being similar in magnitude to that of the FULL experiment. At 14.1 kyr BP, we found that the surface north-south salinity gradient in the FULL experiment is about thirty times the MWF experiment in the SO (Figure 1f). Thus, this increased freshwater flux into the SO area via the SH meltwater forcing enhanced the SO upwelling. Previous studies have also shown that freshwater export from the Antarctic sea ice melt increased the buoyancy of the upwelled SO Circumpolar Deep Water [71]. Consequently, it strengthened the upper branch of the MOC [57,72]. Therefore, this study shows that freshwater discharge to the north and brine rejection to the south of the quasi-permanent Antarctic sea ice boundary steepened the surface north-south salinity gradient, as shown in Figure 1f (signs are reversed). As a result, it strengthened the overturning circulation in the SO.
During the H1 and the YD events, the NH meltwater discharge weakened the AMOC (Figure 3a). As a result, the bipolar seesaw [7][8][9] feedback resulted in deep-ocean warming, increased SH surface air temperature, and Antarctic sea ice retreat [5,78]. However, the AMOC resumed during the ACR, and the bipolar seesaw feedback resulted in the extension of Antarctic sea ice ( Figure 5). In addition, previous experiments have indicated that the rise in austral spring insolation increased the absorption of the incoming shortwave radiation during the last deglaciation. As a result, it warmed the surface-air temperature and further reduced the Antarctic sea ice coverage [75].
The NH summer and the SH spring insolation in the FULL and the ORB experiments increased during the last deglacial period (Figure 1a). The TraCE-21ka FULL and sensitivity experiments indicated that orbital forcing contributed to the deglacial decrease in Antarctic sea ice coverage, as suggested by earlier experiments [98]. However, this study found that the orbital forcing resulted in a relatively milder Antarctic sea ice retreat than the NH meltwater forcing. On the other hand, the NH meltwater forcing showed Antarctic sea ice retreat corresponding to the bipolar seesaw mechanism. Therefore, the TraCE-21ka single-forcing experiments indicated that the combined effect of Northern Hemisphere meltwater flux, increase in insolation and greenhouse gasses, and lowering continental ice sheets could explain the Antarctic sea ice coverage during the last deglacial period.
Recent studies have highlighted that orbital forcing played a substantial role in the changes in AMOC, SO windstress, and SO ventilation during the last deglaciation [76]. Additionally, the authors of [77] provided evidence of orbital-scale changes in the wind-driven SO upwelling. This study shows that the ORB experiment followed the FULL experiment simulated SO upwelling. The ORB experiment also showed a strong correlation between SO upwelling and zonal windstress throughout the last deglacial period. Furthermore, the freshwater discharge from ORB experiment simulated Antarctic sea ice modulated the north-south salinity gradient in the surface of the SO. Thus, this study highlights a significant association of orbital forcing with the evolution of the SO upwelling during the last deglacial period.

Conclusions
The TraCE-21ka ocean model used the eddy parameterization, as it could not resolve eddies. However, we showed that the depth-integrated Eulerian mean upwelling reasonably showed that the deglacial changes in the Southern Ocean upwelling were in agreement with proxy records and numerical simulations. In addition, our results show that the TraCE-21ka simulations agreed with the wind-driven Southern Ocean upwelling hypothesis.
We found that the Northern Hemisphere meltwater discharge contributed to the strength of the upper (AMOC) and lower (Antarctic Bottom Water and its derivatives) branch of the Meridional Overturning Circulation in the Southern Ocean. The meltwater experiment also showed the Antarctic sea ice coverage in response to the bipolar seesaw mechanism compared to a gradual decrease by the orbital experiment throughout the deglacial period. However, the TraCE-21ka orbital experiment, instead of the meltwater experiment, simulated the evolution of the Southern Ocean upwelling reasonably well during the last deglacial period. The orbital experiment also showed a strong correlation between the deglacial Southern Ocean windstress and upwelling, emphasizing the role of wind in Southern Ocean upwelling. Unfortunately, due to data availability, we could not study the seasonal change in the Southern Ocean upwelling. Moreover, it is worth noting that changing continental ice sheets and increasing greenhouse gasses may also contribute to the Southern Ocean dynamics.
Our results show that the coverage of Antarctic quasi-permanent sea ice (Southern Ocean surface area covered with more than eighty percent sea ice) overlapped with the buoyancy loss zone. Thus, a shift in the extent of the deglacial Antarctic sea ice can change the coverage of the surface buoyancy flux. Our results also indicate that Antarctic sea ice and associated surface freshwater discharge played a crucial part in modulating the Southern Ocean upwelling during the last deglaciation.
We found that increased freshwater flux on the surface of the Southern Ocean increased the north-south surface salinity gradient and enhanced the Southern Ocean upwelling. Our study helps understand the physical mechanisms that controlled the Southern Ocean dynamics. The physical mechanisms discussed in our manuscript hold in present-day Southern Ocean dynamics. Presently, we are experiencing ocean warming and Antarctic sea ice melt-induced freshening of the Antarctic shelf waters in response to global warming. Recent studies have also shown that salinity [99,100] and Antarctic meltwater [101] changes in the Southern Ocean. Thus, it is vital to understand the role of freshwater discharge in the Southern Ocean to understand present-day and future climate projections. Our study helps understand the role of freshwater discharge in the Southern Ocean and provides prospects for further scientific development. Thus, our manuscript is vital to understand the physical mechanisms that contribute to the Southern Ocean processes, which can help to understand present-day and future changes. It has to be considered that the impact of changes in the Southern Ocean upwelling cannot precisely explain the variations in atmospheric carbon dioxide concentrations. Besides the Southern Ocean upwelling changes, many other processes and mechanisms, such as ocean biological pump, vertical mixing, properties of the upwelling water, ocean warming, and weathering, can contribute to the deglacial changes in atmospheric carbon dioxide concentrations.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/su14052927/s1: Figure S1: The upper and lower meridional overturning cells during the Heinrich 1 and the onset of the Holocene; Figure S2: The Antarctic sea ice coverage between the Heinrich 1 and the onset of the Holocene;